TRD/qaRec/AliTRDpidRefMakerLQ.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: AliTRDpidRefMakerLQ.cxx 34163 2009-08-07 11:28:51Z cblume $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 //
  21 //  TRD calibration class for building reference data for PID
  22 //  - 2D reference histograms (responsible A.Bercuci)
  23 //  - 3D reference histograms (not yet implemented) (responsible A.Bercuci)
  24 //  - Neural Network (responsible A.Wilk)
  25 //
  26 //   Origin
  27 //   Alex Bercuci  (A.Bercuci@gsi.de)
  28 //
  29 ///////////////////////////////////////////////////////////////////////////////
  30
  31 #include <TFile.h>
  32 #include <TROOT.h>
  33 #include <TTree.h>
  34 #include <TMath.h>
  35 #include <TEventList.h>
  36 #include <TH2D.h>
  37 #include <TH2I.h>
  38 #include <TPrincipal.h>
  39 #include <TVector3.h>
  40 #include <TLinearFitter.h>
  41 #include <TCanvas.h>
  42 #include <TEllipse.h>
  43 #include <TMarker.h>
  44
  45 #include "AliLog.h"
  46 #include "../../STAT/TKDPDF.h"
  47 #include "AliTRDpidRefMakerLQ.h"
  48 #include "../Cal/AliTRDCalPID.h"
  49 #include "AliTRDseedV1.h"
  50 #include "AliTRDcalibDB.h"
  51 #include "AliTRDgeometry.h"
  52
  53 ClassImp(AliTRDpidRefMakerLQ)
  54
  55 //__________________________________________________________________
  56 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
  57   :AliTRDpidRefMaker("PidRefMakerLQ", "PID(LQ) Reference Maker")
  58   ,fPbin(-1)
  59   ,fSbin(-1)
  60 {
  61   //
  62   // AliTRDpidRefMakerLQ default constructor
  63   //
  64
  65   memset(fH2dEdx, 0x0, AliPID::kSPECIES*sizeof(TH2*));
  66 }
  67
  68 //__________________________________________________________________
  69 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
  70 {
  71   //
  72   // AliTRDCalPIDQRef destructor
  73   //
  74
  75   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
  76     if(fH2dEdx[ispec]) delete fH2dEdx[ispec];
  77   }
  78 }
  79
  80 //________________________________________________________________________
  81 void AliTRDpidRefMakerLQ::CreateOutputObjects()
  82 {
  83   // Create histograms
  84   // Called once
  85
  86   AliTRDpidRefMaker::CreateOutputObjects();
  87
  88   fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
  89   fData->Branch("s", &fSbin, "l/b");
  90   fData->Branch("p", &fPbin, "p/b");
  91   fData->Branch("dEdx" , fdEdx , Form("dEdx[%d]/F", GetNslices()));
  92 }
  93
  94
  95 //________________________________________________________________________
  96 Float_t* AliTRDpidRefMakerLQ::GetdEdx(AliTRDseedV1 *trklt)
  97 {
  98   trklt->CookdEdx(AliTRDpidUtil::kLQslices);
  99   const Float_t *dedx = trklt->GetdEdx();
 100   if(dedx[0]+dedx[1] <= 0.) return 0x0;
 101   if(dedx[2] <= 0.) return 0x0;
 102
 103   fdEdx[0] = TMath::Log(dedx[0]+dedx[1]);
 104   fdEdx[1] = TMath::Log(dedx[2]);
 105   return fdEdx;
 106 }
 107
 108
 109 //________________________________________________________________________
 110 void AliTRDpidRefMakerLQ::Fill()
 111 {
 112   // get particle type
 113   fSbin = TMath::LocMax(AliPID::kSPECIES, fPID);
 114   // get momentum bin
 115   fPbin = AliTRDpidUtil::GetMomentumBin(fP);
 116   // fill data
 117   fData->Fill();
 118 }
 119
 120 //________________________________________________________________________
 121 Bool_t AliTRDpidRefMakerLQ::PostProcess()
 122 {
 123   TFile::Open(Form("TRD.Calib%s.root", GetName()));
 124   fData = dynamic_cast<TTree*>(gFile->Get(GetName()));
 125   if (!fData) {
 126     AliError("Tree not available");
 127     return kFALSE;
 128   }
 129   AliDebug(2, Form("Data[%d]", fData->GetEntries()));
 130
 131   Float_t *data[] = {0x0, 0x0};
 132   TPrincipal principal(2, "ND");
 133
 134   TCanvas *cc = new TCanvas("cc", "", 500, 500);
 135
 136   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
 137     for(Int_t is=AliPID::kSPECIES; is--;) {
 138       principal.Clear();
 139       Int_t n = fData->Draw("dEdx[0]:dEdx[1]", Form("p==%d&&s==%d", ip, is), "goff");
 140       AliDebug(2, Form("pBin[%d] sBin[%d] n[%d]", ip, is, n));
 141       if(n<10){
 142         AliWarning(Form("Not enough data for %s[%d].", AliPID::ParticleShortName(is), ip));
 143         continue;
 144       }
 145       // allocate storage
 146       data[0] = new Float_t[n];data[1] = new Float_t[n];
 147
 148       // fill and make principal
 149       Double_t *v1 = fData->GetV1(), *v2 = fData->GetV2();
 150       while(n--){
 151         Double_t dedx[] = {v1[n], v2[n]};
 152         principal.AddRow(dedx);
 153       }
 154       principal.MakePrincipals();
 155
 156 //       // calculate covariance ellipse
 157 //       eValues  = principal.GetEigenValues();
 158 //       x0  = 0.;
 159 //       y0  = 0.;
 160 //       rx  = 3.5*sqrt((*eValues)[0]);
 161 //       ry  = 3.5*sqrt((*eValues)[1]);
 162
 163       // rotate to principal axis
 164       const Double_t *xx; Double_t rxy[2];
 165       while((xx = principal.GetRow(++n))){
 166         principal.X2P(xx, rxy);
 167         data[0][n]=rxy[0]; data[1][n]=rxy[1];
 168         //hProj->Fill(rxy[0], rxy[1]);
 169       }
 170
 171       // estimate acceptable statistical error per bucket
 172
 173       // estimate number of buckets
 174
 175       // build PDF
 176       TKDPDF pdf(n, 2, 10, data);
 177       printf("PDF nodes[%d]\n", pdf.GetNTNodes());
 178       pdf.SetStore();
 179       pdf.SetAlpha(5.);
 180       Float_t *c, v, ve; Double_t r, e;
 181       for(Int_t in=pdf.GetNTNodes(); in--;){
 182         pdf.GetCOGPoint(in, c, v, ve);
 183         rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
 184         printf("%2d x[%f] y[%f]\n", in, rxy[0], rxy[1]);
 185         pdf.Eval(rxy, r, e, kTRUE);
 186       }
 187       pdf.DrawBins(0,1,-6,6,-6,6);cc->Modified(); cc->Update();
 188       AliDebug(2, Form("n[%d]", n));
 189
 190       //pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
 191
 192
 193       delete [] data[0]; delete [] data[1];
 194     }
 195   }
 196
 197   return kTRUE; // testing protection
 198 }
 199
 200
 201 //__________________________________________________________________
 202 void    AliTRDpidRefMakerLQ::Reset()
 203 {
 204   //
 205   // Reset reference histograms
 206   //
 207
 208   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 209     if(fH2dEdx[ispec]) fH2dEdx[ispec]->Reset();
 210   }
 211 }
 212
 213 //__________________________________________________________________
 214 void  AliTRDpidRefMakerLQ::SaveReferences(const Int_t mom, const char *fn)
 215 {
 216   //
 217   // Save the reference histograms
 218   //
 219
 220   TFile *fSave = 0x0;
 221   TListIter it((TList*)gROOT->GetListOfFiles());
 222   Bool_t kFOUND = kFALSE;
 223   TDirectory *pwd = gDirectory;
 224   while((fSave=(TFile*)it.Next()))
 225     if(strcmp(fn, fSave->GetName())==0){
 226       kFOUND = kTRUE;
 227       break;
 228     }
 229   if(!kFOUND) fSave = TFile::Open(fn, "RECREATE");
 230   fSave->cd();
 231
 232   // save dE/dx references
 233   TH2 *h2 = 0x0;
 234   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 235     h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::GetPartSymb(ispec), mom));
 236     h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::GetPartName(ispec), mom));
 237     h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
 238     h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
 239     h2->GetZaxis()->SetTitle("Entries");
 240     h2->Write();
 241   }
 242
 243   pwd->cd();
 244 }
 245