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 //Root includes
  19 #include <TObjArray.h>
  20 #include <TH2F.h>
  21 #include <TString.h>
  22 #include <TVectorD.h>
  23 #include <TAxis.h>
  24 #include <TLinearFitter.h>
  25 #include <TMath.h>
  26 #include <TTreeStream.h>
  27
  28 //header file
  29 #include "AliTRDCalibraVdriftLinearFit.h"
  30
  31 ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
  32
  33 //_____________________________________________________________________
  34 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
  35   TObject(),
  36   fVersion(0),
  37   fLinearFitterHistoArray(540),
  38   fLinearFitterPArray(540),
  39   fLinearFitterEArray(540)
  40 {
  41     //
  42     // default constructor
  43     //
  44 }
  45 //_____________________________________________________________________
  46 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
  47   TObject(ped),
  48   fVersion(ped.fVersion),
  49   fLinearFitterHistoArray(540),
  50   fLinearFitterPArray(540),
  51   fLinearFitterEArray(540)
  52 {
  53     //
  54     // copy constructor
  55     //
  56   for (Int_t idet = 0; idet < 540; idet++){
  57     const TVectorD     *vectorE     = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
  58     const TVectorD     *vectorP     = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
  59     const TH2F         *hped        = (TH2F*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
  60
  61     if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
  62     if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
  63     if ( hped != 0x0 ){
  64       TH2F *hNew = new TH2F(*hped);
  65       hNew->SetDirectory(0);
  66       fLinearFitterHistoArray.AddAt(hNew,idet);
  67     }
  68
  69   }
  70 }
  71 //_____________________________________________________________________
  72 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
  73   TObject(),
  74   fVersion(0),
  75   fLinearFitterHistoArray(540),
  76   fLinearFitterPArray(540),
  77   fLinearFitterEArray(540)
  78 {
  79   //
  80   // constructor from a TObjArray
  81   //
  82   for (Int_t idet = 0; idet < 540; idet++){
  83     const TH2F         *hped        = (TH2F*)obja.UncheckedAt(idet);
  84     if ( hped != 0x0 ){
  85       TH2F *hNew = new TH2F(*hped);
  86       hNew->SetDirectory(0);
  87       fLinearFitterHistoArray.AddAt(hNew,idet);
  88     }
  89   }
  90 }
  91 //_____________________________________________________________________
  92 AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const  AliTRDCalibraVdriftLinearFit &source)
  93 {
  94   //
  95   // assignment operator
  96   //
  97   if (&source == this) return *this;
  98   new (this) AliTRDCalibraVdriftLinearFit(source);
  99
 100   return *this;
 101 }
 102 //_____________________________________________________________________
 103 AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
 104 {
 105   //
 106   // destructor
 107   //
 108 }
 109 //_____________________________________________________________________
 110 void AliTRDCalibraVdriftLinearFit::Add(AliTRDCalibraVdriftLinearFit *ped)
 111 {
 112   //
 113   // Add histo
 114   //
 115
 116   fVersion++;
 117
 118   for (Int_t idet = 0; idet < 540; idet++){
 119     const TH2F         *hped        = (TH2F*)ped->GetLinearFitterHisto(idet);
 120     //printf("idet %d\n",idet);
 121     if ( hped != 0x0 ){
 122       //printf("add\n");
 123       TH2F *hped1 = GetLinearFitterHisto(idet,kTRUE);
 124       //printf("test2\n");
 125       hped1->SetDirectory(0);
 126       //printf("test4\n");
 127       hped1->Add(hped);
 128       //printf("test3\n");
 129     }
 130   }
 131 }
 132 //______________________________________________________________________________________
 133 TH2F* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
 134 {
 135     //
 136     // return pointer to TH2F histo
 137     // if force is true create a new histo if it doesn't exist allready
 138     //
 139     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
 140         return (TH2F*)fLinearFitterHistoArray.UncheckedAt(detector);
 141
 142     // if we are forced and TLinearFitter doesn't yes exist create it
 143
 144     // new TH2F
 145     TString name("LFDV");
 146     name += detector;
 147     name += "version";
 148     name +=  fVersion;
 149
 150     TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name
 151                           ,100,-1.0,1.0,100
 152                           ,-2.0,2.0);
 153     lfdv->SetXTitle("tan(phi_{track})");
 154     lfdv->SetYTitle("dy/dt");
 155     lfdv->SetZTitle("Number of clusters");
 156     lfdv->SetStats(0);
 157     lfdv->SetDirectory(0);
 158
 159     fLinearFitterHistoArray.AddAt(lfdv,detector);
 160     return lfdv;
 161 }
 162 //______________________________________________________________________________________
 163 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
 164 {
 165     //
 166     // return param for this detector
 167     //
 168   if ( fLinearFitterPArray.UncheckedAt(detector) ){
 169     const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
 170     if(!param) param = new TVectorD(2);
 171     for(Int_t k = 0; k < 2; k++){
 172       (*param)[k] = (*vectorP)[k];
 173     }
 174     return kTRUE;
 175   }
 176   else return kFALSE;
 177
 178 }
 179 //______________________________________________________________________________________
 180 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
 181 {
 182     //
 183     // return error for this detector
 184     //
 185   if ( fLinearFitterEArray.UncheckedAt(detector) ){
 186     const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
 187     if(!error) error = new TVectorD(3);
 188     for(Int_t k = 0; k < 3; k++){
 189       (*error)[k] = (*vectorE)[k];
 190     }
 191     return kTRUE;
 192   }
 193   else return kFALSE;
 194
 195 }
 196 //______________________________________________________________________________________
 197 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
 198 {
 199     //
 200     // Fill the 2D histos for debugging
 201     //
 202
 203   ((TH2F *) GetLinearFitterHisto(detector,kTRUE))->Fill(tnp,pars1);
 204
 205 }
 206 //____________Functions fit Online CH2d________________________________________
 207 void AliTRDCalibraVdriftLinearFit::FillPEArray()
 208 {
 209   //
 210   // Fill fLinearFitterPArray and fLinearFitterEArray from inside
 211   //
 212
 213
 214   Int_t *arrayI = new Int_t[540];
 215   for(Int_t k = 0; k< 540; k++){
 216     arrayI[k] = 0;
 217   }
 218
 219   // Loop over histos
 220   for(Int_t cb = 0; cb < 540; cb++){
 221     const TH2F *linearfitterhisto = (TH2F*)fLinearFitterHistoArray.UncheckedAt(cb);
 222     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
 223
 224     if ( linearfitterhisto != 0 ){
 225
 226       // Fill a linearfitter
 227       TAxis *xaxis = linearfitterhisto->GetXaxis();
 228       TAxis *yaxis = linearfitterhisto->GetYaxis();
 229       TLinearFitter linearfitter = TLinearFitter(2,"pol1");
 230       //printf("test\n");
 231       for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
 232         for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
 233           if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
 234             Double_t x = xaxis->GetBinCenter(ibinx+1);
 235             Double_t y = yaxis->GetBinCenter(ibiny+1);
 236             for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
 237               linearfitter.AddPoint(&x,y);
 238               arrayI[cb]++;
 239             }
 240           }
 241         }
 242       }
 243
 244       //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
 245
 246       // Eval the linear fitter
 247       if(arrayI[cb]>10){
 248         TVectorD  *par  = new TVectorD(2);
 249         TVectorD   pare = TVectorD(2);
 250         TVectorD  *parE = new TVectorD(3);
 251         linearfitter.Eval();
 252         linearfitter.GetParameters(*par);
 253         linearfitter.GetErrors(pare);
 254         Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
 255         (*parE)[0] = pare[0]*ppointError;
 256         (*parE)[1] = pare[1]*ppointError;
 257         (*parE)[2] = (Double_t) arrayI[cb];
 258         fLinearFitterPArray.AddAt(par,cb);
 259         fLinearFitterEArray.AddAt(parE,cb);
 260         //par->Print();
 261         //parE->Print();
 262       }
 263     }// if something
 264   }
 265
 266 }