EMCAL/AliEMCALPID.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 /* History of cvs commits:
  18  *
  19  * $Log$
  20  *
  21  */
  22 //    to compute PID for all the clusters in ESDs.root file
  23 //     the ESDs.root have to be in the same directory as the class
  24 //
  25 //
  26 //
  27 //
  28 //
  29 //      AliEMCALPID::CalculPID(Energy,Lambda0)
  30 //        Calcul PID for all clusters in AliESDs.root file
  31 //        keep this function for the moment for a simple verification, could be removed
  32 //
  33 //
  34 //
  35 //   AliEMCALPID::CalculPID(Energy,Lambda0)
  36 //    calcul PID Weght for a cluster with Energy, Lambda0 .
  37 //    Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
  38 //
  39 //
  40 //
  41 //      kElectron :  fPIDFinal[0]
  42 //      kMuon     :  fPIDFinal[1]
  43 //      kPion     :  fPIDFinal[2]
  44 //      kKaon     :  fPIDFinal[3]
  45 //      kProton   :  fPIDFinal[4]
  46 //      kPhoton   :  fPIDFinal[5]
  47 //      kPi0      :  fPIDFinal[6]
  48 //      kNeutron  :  fPIDFinal[7]
  49 //      kKaon0    :  fPIDFinal[8]
  50 //      kEleCon   :  fPIDFinal[9]
  51 //      kUnknown  :  fPIDFinal[10]
  52 //
  53 //
  54 //    PID[3] is a simple PID for
  55 //      Electron & Photon  PID[0]
  56 //                    Pi0  PID[1]
  57 //                 Hadron  PID[2]
  58 //
  59 //
  60 //
  61 //
  62 //
  63 // --- ROOT system ---
  64
  65 // standard C++ includes
  66 #include <Riostream.h>
  67 // #include <cstdlib>
  68 // #include <cmath>
  69
  70 // ROOT includes
  71 #include "TTree.h"
  72 #include "TStyle.h"
  73 #include "TVector3.h"
  74 #include "TBranch.h"
  75 #include "TClonesArray.h"
  76 #include "TCanvas.h"
  77 #include "TLorentzVector.h"
  78 #include "TMath.h"
  79 #include "TFile.h"
  80 #include "TH1.h"
  81 #include "TH2.h"
  82 #include "TParticle.h"
  83
  84 // // STEER includes
  85 // #include "AliRun.h"
  86 // #include "AliRunLoader.h"
  87 // #include "AliHeader.h"
  88 // #include "AliLoader.h"
  89 // #include "AliStack.h"
  90 // #include "AliESDtrack.h"
  91 // #include "AliESD.h"
  92 #include "AliLog.h"
  93 #include "AliEMCALPID.h"
  94
  95 ClassImp(AliEMCALPID)
  96
  97 AliEMCALPID::AliEMCALPID()
  98 {
  99 //
 100 // Constructor.
 101 // Initialize all constant values which have to be used
 102 // during PID algorithm execution
 103 //
 104
 105         // set flag for printing to FALSE by default
 106         fPrintInfo = kFALSE;
 107
 108         // as a first step, all array elements are initialized to 0.0
 109         Int_t i, j;
 110         for (i = 0; i < 6; i++) {
 111                 for (j = 0; j < 6; j++) {
 112                         fGamma[i][j] = fHadron[i][j] = fPiZero5to10[i][j] = fPiZero10to60[i][j] = 0.;
 113                 }
 114         }
 115
 116         // then, only the ones which must be not zero are initialized
 117         // while the others will remain to the value 0.0
 118
 119         fGamma[0][0] =  0.038022;
 120         fGamma[0][1] = -0.0001883;
 121         fGamma[0][2] =  5.449e-06;
 122
 123         fGamma[1][0] =  0.207313;
 124         fGamma[1][1] = -0.000978;
 125         fGamma[1][2] =  0.00001634;
 126
 127         fGamma[2][0] =  0.043364;
 128         fGamma[2][1] = -0.0002048;
 129         fGamma[2][2] =  8.661e-06;
 130         fGamma[2][3] = -1.353e-07;
 131
 132         fGamma[3][0] =  0.265004;
 133         fGamma[3][1] =  0.061298;
 134         fGamma[3][2] = -0.003203;
 135         fGamma[3][3] =  4.73e-05;
 136
 137         fGamma[4][0] =  0.243579;
 138         fGamma[4][1] = -1.614e-05;
 139
 140         fGamma[5][0] =  0.002942;
 141         fGamma[5][1] = -3.976e-05;
 142
 143         fHadron[0][0] =  0.011945 / 3.;
 144         fHadron[0][1] =  0.000386 / 3.;
 145         fHadron[0][2] = -0.000014 / 3.;
 146         fHadron[0][3] =  1.336e-07 / 3.;
 147
 148         fHadron[1][0] =  0.496544;
 149         fHadron[1][1] = -0.003226;
 150         fHadron[1][2] =  0.00001678;
 151
 152         fHadron[2][0] =  0.144838;
 153         fHadron[2][1] = -0.002954;
 154         fHadron[2][2] =  0.00008754;
 155         fHadron[2][3] = -7.587e-07;
 156
 157         fHadron[3][0] =  1.264461 / 7.;
 158         fHadron[3][1] =  0.002097 / 7.;
 159
 160         fHadron[4][0] =  0.261950;
 161         fHadron[4][1] = -0.001078;
 162         fHadron[4][2] =  0.00003237;
 163         fHadron[4][3] = -3.241e-07;
 164         fHadron[4][4] =  0.;
 165         fHadron[4][5] =  0.;
 166         fHadron[5][0] =  0.010317;
 167         fHadron[5][1] =  0.;
 168         fHadron[5][2] =  0.;
 169         fHadron[5][3] =  0.;
 170         fHadron[5][4] =  0.;
 171         fHadron[5][5] =  0.;
 172
 173         fPiZero5to10[0][0] = 0.009138;
 174         fPiZero5to10[0][1] = 0.0006377;
 175
 176         fPiZero5to10[1][0] = 0.08;
 177
 178         fPiZero5to10[2][0] = -0.061119;
 179         fPiZero5to10[2][1] =  0.019013;
 180
 181         fPiZero5to10[3][0] =  0.2;
 182
 183         fPiZero5to10[4][0] =  0.252044;
 184         fPiZero5to10[4][1] = -0.002315;
 185
 186         fPiZero5to10[5][0] =  0.002942;
 187         fPiZero5to10[5][1] = -3.976e-05;
 188
 189         fPiZero10to60[0][0] =  0.009138;
 190         fPiZero10to60[0][1] =  0.0006377;
 191
 192         fPiZero10to60[1][0] =  1.272837;
 193         fPiZero10to60[1][1] = -0.069708;
 194         fPiZero10to60[1][2] =  0.001568;
 195         fPiZero10to60[1][3] = -1.162e-05;
 196
 197         fPiZero10to60[2][0] =  0.139703;
 198         fPiZero10to60[2][1] =  0.003687;
 199         fPiZero10to60[2][2] = -0.000568;
 200         fPiZero10to60[2][3] =  1.498e-05;
 201         fPiZero10to60[2][4] = -1.174e-07;
 202
 203         fPiZero10to60[3][0] = -0.826367;
 204         fPiZero10to60[3][1] =  0.096951;
 205         fPiZero10to60[3][2] = -0.002215;
 206         fPiZero10to60[3][3] =  2.523e-05;
 207
 208         fPiZero10to60[4][0] =  0.249890;
 209         fPiZero10to60[4][1] = -0.000063;
 210
 211         fPiZero10to60[5][0] =  0.002942;
 212         fPiZero10to60[5][1] = -3.976e-05;
 213
 214         fPIDWeight[0] = -1;
 215         fPIDWeight[1] = -1;
 216         fPIDWeight[2] = -1;
 217         fReconstructor = kFALSE;
 218 }
 219 //
 220 //
 221 void AliEMCALPID::RunPID(AliESD *esd)
 222 {
 223 //
 224 // Make the PID for all the EMCAL clusters containedin the ESDs File
 225 // but just gamma/PiO/Hadron
 226 //
 227         // trivial check against NULL object passed
 228
 229   if (esd == 0x0) {
 230     AliInfo("NULL ESD object passed!!" );
 231     return ;
 232   }
 233   Int_t nClusters = esd->GetNumberOfEMCALClusters();
 234   Int_t firstCluster = esd->GetFirstEMCALCluster();
 235   Double_t energy, lambda0;
 236   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
 237
 238     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
 239     energy = clust->GetClusterEnergy();
 240     lambda0 = clust->GetM02();
 241     // verify cluster type
 242     Int_t clusterType= clust->GetClusterType();
 243     if (clusterType == AliESDCaloCluster::kClusterv1 && lambda0 != 0 && energy > 5 && energy < 1000) {
 244       // reject clusters with lambda0 = 0
 245       // reject clusters with energy < 5 GeV
 246       ComputePID(energy, lambda0);
 247       if (fPrintInfo) {
 248         AliInfo("___________________________________________________");
 249         AliInfo(Form( "Particle Energy = %f",energy));
 250         AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
 251         AliInfo("PIDWeight of the particle :" );
 252         AliInfo(Form( " GAMMA  : %f",fPID[0] ));
 253         AliInfo(Form( " PiZero : %f",fPID[1] ));
 254         AliInfo(Form( " HADRON : %f", fPID[2] ));
 255         AliInfo("_________________________________________");
 256         AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
 257         AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
 258         AliInfo(Form( " kPion       : %f", fPIDFinal[2] ));
 259         AliInfo(Form( " kKaon       : %f", fPIDFinal[3] ));
 260         AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
 261         AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
 262         AliInfo(Form( " kPi0        : %f", fPIDFinal[6] ));
 263         AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
 264         AliInfo(Form( " kKaon0  : %f", fPIDFinal[8] ));
 265         AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
 266         AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
 267         AliInfo("___________________________________________________");
 268       }
 269       if(fReconstructor) // In case it is called during reconstruction.
 270         clust->SetPid(fPIDFinal);
 271     } // end if (clusterType...)
 272   } // end for (iCluster...)
 273 }
 274 //
 275 //
 276 void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
 277 {
 278 //
 279 // This is the main command, which uses the distributions computed and parametrised,
 280 // and gives the PID by the bayesian method.
 281 //
 282         TArrayD paramDistribGamma  = DistLambda0(energy, 1);
 283         TArrayD paramDistribPiZero = DistLambda0(energy, 2);
 284         TArrayD paramDistribHadron = DistLambda0(energy, 3);
 285
 286         Bool_t norm = kFALSE;
 287
 288         fProbGamma   = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
 289         fProbGamma  += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
 290         fProbPiZero  = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
 291         fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
 292         fProbHadron  = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
 293         fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
 294
 295         // compute PID Weight
 296         fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
 297         fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
 298         fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
 299
 300         SetPID(fPIDWeight[0], 0);
 301         SetPID(fPIDWeight[1], 1);
 302         SetPID(fPIDWeight[2], 2);
 303
 304         // sortie ecran pid Weight only for control (= in english ???)
 305         if (fPrintInfo) {
 306           AliInfo(Form( "Energy in loop = %f", energy) );
 307           AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
 308           AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
 309                 //      AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
 310           AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
 311           AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
 312           AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f",  fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
 313           AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
 314           AliInfo("********************************************************" );
 315         }
 316
 317         fPIDFinal[0]  = fPIDWeight[0]/2;
 318         fPIDFinal[1]  = fPIDWeight[2]/8;
 319         fPIDFinal[2]  = fPIDWeight[2]/8;
 320         fPIDFinal[3]  = fPIDWeight[2]/8;
 321         fPIDFinal[4]  = fPIDWeight[2]/8;
 322         fPIDFinal[5]  = fPIDWeight[0]/2;
 323         fPIDFinal[6]  = fPIDWeight[1]/2;
 324         fPIDFinal[7]  = fPIDWeight[2]/8;
 325         fPIDFinal[8]  = fPIDWeight[2]/8;
 326         fPIDFinal[9]  = fPIDWeight[2]/8;
 327         fPIDFinal[10] = fPIDWeight[2]/8;
 328         fPIDFinal[11] = 0;
 329 }
 330 //
 331 //
 332 TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
 333 {
 334 //
 335 // Compute the values of the parametrised distributions using the data initialised before.
 336 //
 337         Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
 338         Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
 339         TArrayD  distributionParam(6);
 340
 341         switch (type) {
 342                 case 1:
 343                         constGauss  = Polynomial(energy, fGamma[0]);
 344                         meanGauss   = Polynomial(energy, fGamma[1]);
 345                         sigmaGauss  = Polynomial(energy, fGamma[2]);
 346                         constLandau = Polynomial(energy, fGamma[3]);
 347                         mpvLandau   = Polynomial(energy, fGamma[4]);
 348                         sigmaLandau = Polynomial(energy, fGamma[5]);
 349                         break;
 350                 case 2:
 351                         if (energy < 10) {
 352                                 constGauss  = Polynomial(energy, fPiZero5to10[0]);
 353                                 meanGauss   = Polynomial(energy, fPiZero5to10[1]);
 354                                 sigmaGauss  = Polynomial(energy, fPiZero5to10[2]);
 355                                 constLandau = Polynomial(energy, fPiZero5to10[3]);
 356                                 mpvLandau   = Polynomial(energy, fPiZero5to10[4]);
 357                                 sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
 358                         }
 359                         else {
 360                                 constGauss  = Polynomial(energy, fPiZero10to60[0]);
 361                                 meanGauss   = Polynomial(energy, fPiZero10to60[1]);
 362                                 sigmaGauss  = Polynomial(energy, fPiZero10to60[2]);
 363                                 constLandau = Polynomial(energy, fPiZero10to60[3]);
 364                                 mpvLandau   = Polynomial(energy, fPiZero10to60[4]);
 365                                 sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
 366                         }
 367                         break;
 368                 case 3:
 369                         constGauss  = Polynomial(energy, fHadron[0]);
 370                         meanGauss   = Polynomial(energy, fHadron[1]);
 371                         sigmaGauss  = Polynomial(energy, fHadron[2]);
 372                         constLandau = Polynomial(energy, fHadron[3]);
 373                         mpvLandau   = Polynomial(energy, fHadron[4]);
 374                         sigmaLandau = Polynomial(energy, fHadron[5]);
 375                         break;
 376         }
 377
 378         distributionParam[0] = constGauss;
 379         distributionParam[1] = meanGauss;
 380         distributionParam[2] = sigmaGauss;
 381         distributionParam[3] = constLandau;
 382         distributionParam[4] = mpvLandau;
 383         distributionParam[5] = sigmaLandau;
 384
 385         return distributionParam;
 386 }
 387 //
 388 //
 389 Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
 390 {
 391 //
 392 // Compute a polynomial for a given value of 'x'
 393 // with the array of parameters passed as the second arg
 394 //
 395         Double_t y;
 396         y  = params[0];
 397         y += params[1] * x;
 398         y += params[2] * x * x;
 399         y += params[3] * x * x * x;
 400         y += params[4] * x * x * x * x;
 401         y += params[5] * x * x * x * x * x;
 402
 403         return y;
 404 }