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