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