| | 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 | /* History of cvs commits: |
| | 19 | * |
| | 20 | * $Log$ |
| | 21 | * Revision 1.101 2005/05/28 14:19:04 schutz |
| | 22 | * Compilation warnings fixed by T.P. |
| | 23 | * |
| | 24 | */ |
| | 25 | |
| | 26 | //_________________________________________________________________________ |
| | 27 | // Implementation version v1 of the PHOS particle identifier |
| | 28 | // Particle identification based on the |
| | 29 | // - RCPV: distance from CPV recpoint to EMCA recpoint. |
| | 30 | // - TOF |
| | 31 | // - PCA: Principal Components Analysis.. |
| | 32 | // The identified particle has an identification number corresponding |
| | 33 | // to a 9 bits number: |
| | 34 | // -Bit 0 to 2: bit set if RCPV > CpvEmcDistance (each bit corresponds |
| | 35 | // to a different efficiency-purity point of the photon identification) |
| | 36 | // -Bit 3 to 5: bit set if TOF < TimeGate (each bit corresponds |
| | 37 | // to a different efficiency-purity point of the photon identification) |
| | 38 | // -Bit 6 to 9: bit set if Principal Components are |
| | 39 | // inside an ellipse defined by the parameters a, b, c, x0 and y0. |
| | 40 | // (each bit corresponds to a different efficiency-purity point of the |
| | 41 | // photon identification) |
| | 42 | // The PCA (Principal components analysis) needs a file that contains |
| | 43 | // a previous analysis of the correlations between the particles. This |
| | 44 | // file is $ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root. Analysis done for |
| | 45 | // energies between 0.5 and 100 GeV. |
| | 46 | // A calibrated energy is calculated. The energy of the reconstructed |
| | 47 | // cluster is corrected with the formula A + B * E + C * E^2, whose |
| | 48 | // parameters where obtained through the study of the reconstructed |
| | 49 | // energy distribution of monoenergetic photons. |
| | 50 | // |
| | 51 | // All the parameters (RCPV(2 rows-3 columns),TOF(1r-3c),PCA(5r-4c) |
| | 52 | // and calibration(1r-3c))are stored in a file called |
| | 53 | // $ALICE_ROOT/PHOS/Parameters.dat. Each time that AliPHOSPIDv1 is |
| | 54 | // initialized, this parameters are copied to a Matrix (9,4), a |
| | 55 | // TMatrixD object. |
| | 56 | // |
| | 57 | // use case: |
| | 58 | // root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root") |
| | 59 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
| | 60 | // // reading headers from file galice1.root and create RecParticles |
| | 61 | // TrackSegments and RecPoints are used |
| | 62 | // // set file name for the branch RecParticles |
| | 63 | // root [1] p->ExecuteTask("deb all time") |
| | 64 | // // available options |
| | 65 | // // "deb" - prints # of reconstructed particles |
| | 66 | // // "deb all" - prints # and list of RecParticles |
| | 67 | // // "time" - prints benchmarking results |
| | 68 | // |
| | 69 | // root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1",kTRUE) |
| | 70 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
| | 71 | // //Split mode. |
| | 72 | // root [3] p2->ExecuteTask() |
| | 73 | // |
| | 74 | |
| | 75 | |
| | 76 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & |
| | 77 | // Gustavo Conesa April 2002 |
| | 78 | // PCA redesigned by Gustavo Conesa October 2002: |
| | 79 | // The way of using the PCA has changed. Instead of 2 |
| | 80 | // files with the PCA, each one with different energy ranges |
| | 81 | // of application, we use the wide one (0.5-100 GeV), and instead |
| | 82 | // of fixing 3 ellipses for different ranges of energy, it has been |
| | 83 | // studied the dependency of the ellipses parameters with the |
| | 84 | // energy, and they are implemented in the code as a funtion |
| | 85 | // of the energy. |
| | 86 | // |
| | 87 | // |
| | 88 | // |
| | 89 | // --- ROOT system --- |
| | 90 | |
| | 91 | |
| | 92 | // --- Standard library --- |
| | 93 | #include <TMatrixF.h> |
| | 94 | #include "TFormula.h" |
| | 95 | #include "TBenchmark.h" |
| | 96 | #include "TPrincipal.h" |
| | 97 | #include "TFile.h" |
| | 98 | #include "TSystem.h" |
| | 99 | |
| | 100 | // --- AliRoot header files --- |
| | 101 | //#include "AliLog.h" |
| | 102 | #include "AliGenerator.h" |
| | 103 | #include "AliPHOS.h" |
| | 104 | #include "AliPHOSPIDv1.h" |
| | 105 | #include "AliPHOSGetter.h" |
| | 106 | |
| | 107 | ClassImp( AliPHOSPIDv1) |
| | 108 | |
| | 109 | //____________________________________________________________________________ |
| | 110 | AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() |
| | 111 | { |
| | 112 | // default ctor |
| | 113 | |
| | 114 | InitParameters() ; |
| | 115 | fDefaultInit = kTRUE ; |
| | 116 | } |
| | 117 | |
| | 118 | //____________________________________________________________________________ |
| | 119 | AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ):AliPHOSPID(pid) |
| | 120 | { |
| | 121 | // ctor |
| | 122 | InitParameters() ; |
| | 123 | Init() ; |
| | 124 | |
| | 125 | } |
| | 126 | |
| | 127 | //____________________________________________________________________________ |
| | 128 | AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName):AliPHOSPID(alirunFileName, eventFolderName) |
| | 129 | { |
| | 130 | //ctor with the indication on where to look for the track segments |
| | 131 | |
| | 132 | InitParameters() ; |
| | 133 | Init() ; |
| | 134 | fDefaultInit = kFALSE ; |
| | 135 | } |
| | 136 | |
| | 137 | //____________________________________________________________________________ |
| | 138 | AliPHOSPIDv1::~AliPHOSPIDv1() |
| | 139 | { |
| | 140 | // dtor |
| | 141 | fPrincipalPhoton = 0; |
| | 142 | fPrincipalPi0 = 0; |
| | 143 | |
| | 144 | delete [] fX ; // Principal input |
| | 145 | delete [] fPPhoton ; // Photon Principal components |
| | 146 | delete [] fPPi0 ; // Pi0 Principal components |
| | 147 | |
| | 148 | delete fParameters; |
| | 149 | delete fTFphoton; |
| | 150 | delete fTFpiong; |
| | 151 | delete fTFkaong; |
| | 152 | delete fTFkaonl; |
| | 153 | delete fTFhhadrong; |
| | 154 | delete fTFhhadronl; |
| | 155 | delete fDFmuon; |
| | 156 | } |
| | 157 | //____________________________________________________________________________ |
| | 158 | const TString AliPHOSPIDv1::BranchName() const |
| | 159 | { |
| | 160 | |
| | 161 | return GetName() ; |
| | 162 | } |
| | 163 | |
| | 164 | //____________________________________________________________________________ |
| | 165 | void AliPHOSPIDv1::Init() |
| | 166 | { |
| | 167 | // Make all memory allocations that are not possible in default constructor |
| | 168 | // Add the PID task to the list of PHOS tasks |
| | 169 | |
| | 170 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
| | 171 | if(!gime) |
| | 172 | gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ; |
| | 173 | |
| | 174 | if ( !gime->PID() ) |
| | 175 | gime->PostPID(this) ; |
| | 176 | } |
| | 177 | |
| | 178 | //____________________________________________________________________________ |
| | 179 | void AliPHOSPIDv1::InitParameters() |
| | 180 | { |
| | 181 | // Initialize PID parameters |
| | 182 | fWrite = kTRUE ; |
| | 183 | fRecParticlesInRun = 0 ; |
| | 184 | fNEvent = 0 ; |
| | 185 | fRecParticlesInRun = 0 ; |
| | 186 | fBayesian = kTRUE ; |
| | 187 | SetParameters() ; // fill the parameters matrix from parameters file |
| | 188 | SetEventRange(0,-1) ; |
| | 189 | |
| | 190 | // initialisation of response function parameters |
| | 191 | // Tof |
| | 192 | |
| | 193 | // // Photons |
| | 194 | // fTphoton[0] = 0.218 ; |
| | 195 | // fTphoton[1] = 1.55E-8 ; |
| | 196 | // fTphoton[2] = 5.05E-10 ; |
| | 197 | // fTFphoton = new TFormula("ToF response to photons" , "gaus") ; |
| | 198 | // fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; |
| | 199 | |
| | 200 | // // Pions |
| | 201 | // //Gaus (0 to max probability) |
| | 202 | // fTpiong[0] = 0.0971 ; |
| | 203 | // fTpiong[1] = 1.58E-8 ; |
| | 204 | // fTpiong[2] = 5.69E-10 ; |
| | 205 | // fTFpiong = new TFormula("ToF response to pions" , "gaus") ; |
| | 206 | // fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; |
| | 207 | |
| | 208 | // // Kaons |
| | 209 | // //Gaus (0 to max probability) |
| | 210 | // fTkaong[0] = 0.0542 ; |
| | 211 | // fTkaong[1] = 1.64E-8 ; |
| | 212 | // fTkaong[2] = 6.07E-10 ; |
| | 213 | // fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; |
| | 214 | // fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; |
| | 215 | // //Landau (max probability to inf) |
| | 216 | // fTkaonl[0] = 0.264 ; |
| | 217 | // fTkaonl[1] = 1.68E-8 ; |
| | 218 | // fTkaonl[2] = 4.10E-10 ; |
| | 219 | // fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; |
| | 220 | // fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; |
| | 221 | |
| | 222 | // //Heavy Hadrons |
| | 223 | // //Gaus (0 to max probability) |
| | 224 | // fThhadrong[0] = 0.0302 ; |
| | 225 | // fThhadrong[1] = 1.73E-8 ; |
| | 226 | // fThhadrong[2] = 9.52E-10 ; |
| | 227 | // fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; |
| | 228 | // fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; |
| | 229 | // //Landau (max probability to inf) |
| | 230 | // fThhadronl[0] = 0.139 ; |
| | 231 | // fThhadronl[1] = 1.745E-8 ; |
| | 232 | // fThhadronl[2] = 1.00E-9 ; |
| | 233 | // fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; |
| | 234 | // fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; |
| | 235 | |
| | 236 | // Photons |
| | 237 | fTphoton[0] = 7.83E8 ; |
| | 238 | fTphoton[1] = 1.55E-8 ; |
| | 239 | fTphoton[2] = 5.09E-10 ; |
| | 240 | fTFphoton = new TFormula("ToF response to photons" , "gaus") ; |
| | 241 | fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; |
| | 242 | |
| | 243 | // Pions |
| | 244 | //Gaus (0 to max probability) |
| | 245 | fTpiong[0] = 6.73E8 ; |
| | 246 | fTpiong[1] = 1.58E-8 ; |
| | 247 | fTpiong[2] = 5.87E-10 ; |
| | 248 | fTFpiong = new TFormula("ToF response to pions" , "gaus") ; |
| | 249 | fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; |
| | 250 | |
| | 251 | // Kaons |
| | 252 | //Gaus (0 to max probability) |
| | 253 | fTkaong[0] = 3.93E8 ; |
| | 254 | fTkaong[1] = 1.64E-8 ; |
| | 255 | fTkaong[2] = 6.07E-10 ; |
| | 256 | fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; |
| | 257 | fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; |
| | 258 | //Landau (max probability to inf) |
| | 259 | fTkaonl[0] = 2.0E9 ; |
| | 260 | fTkaonl[1] = 1.68E-8 ; |
| | 261 | fTkaonl[2] = 4.10E-10 ; |
| | 262 | fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; |
| | 263 | fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; |
| | 264 | |
| | 265 | //Heavy Hadrons |
| | 266 | //Gaus (0 to max probability) |
| | 267 | fThhadrong[0] = 2.02E8 ; |
| | 268 | fThhadrong[1] = 1.73E-8 ; |
| | 269 | fThhadrong[2] = 9.52E-10 ; |
| | 270 | fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; |
| | 271 | fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; |
| | 272 | //Landau (max probability to inf) |
| | 273 | fThhadronl[0] = 1.10E9 ; |
| | 274 | fThhadronl[1] = 1.74E-8 ; |
| | 275 | fThhadronl[2] = 1.00E-9 ; |
| | 276 | fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; |
| | 277 | fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; |
| | 278 | |
| | 279 | |
| | 280 | |
| | 281 | // Shower shape: dispersion gaussian parameters |
| | 282 | // Photons |
| | 283 | |
| | 284 | // fDphoton[0] = 4.62e-2; fDphoton[1] = 1.39e-2 ; fDphoton[2] = -3.80e-2;//constant |
| | 285 | // fDphoton[3] = 1.53 ; fDphoton[4] =-6.62e-2 ; fDphoton[5] = 0.339 ;//mean |
| | 286 | // fDphoton[6] = 6.89e-2; fDphoton[7] =-6.59e-2 ; fDphoton[8] = 0.194 ;//sigma |
| | 287 | |
| | 288 | // fDpi0[0] = 0.0586 ; fDpi0[1] = 1.06E-3 ; fDpi0[2] = 0. ;//constant |
| | 289 | // fDpi0[3] = 2.67 ; fDpi0[4] =-2.00E-2 ; fDpi0[5] = 9.37E-5 ;//mean |
| | 290 | // fDpi0[6] = 0.153 ; fDpi0[7] = 9.34E-4 ; fDpi0[8] =-1.49E-5 ;//sigma |
| | 291 | |
| | 292 | // fDhadron[0] = 1.61E-2 ; fDhadron[1] = 3.03E-3 ; fDhadron[2] = 1.01E-2 ;//constant |
| | 293 | // fDhadron[3] = 3.81 ; fDhadron[4] = 0.232 ; fDhadron[5] =-1.25 ;//mean |
| | 294 | // fDhadron[6] = 0.897 ; fDhadron[7] = 0.0987 ; fDhadron[8] =-0.534 ;//sigma |
| | 295 | |
| | 296 | fDphoton[0] = 1.5 ; fDphoton[1] = 0.49 ; fDphoton[2] =-1.7E-2 ;//constant |
| | 297 | fDphoton[3] = 1.5 ; fDphoton[4] = 4.0E-2 ; fDphoton[5] = 0.21 ;//mean |
| | 298 | fDphoton[6] = 4.8E-2 ; fDphoton[7] =-0.12 ; fDphoton[8] = 0.27 ;//sigma |
| | 299 | fDphoton[9] = 16.; //for E> fDphoton[9] parameters calculated at fDphoton[9] |
| | 300 | |
| | 301 | fDpi0[0] = 0.25 ; fDpi0[1] = 3.3E-2 ; fDpi0[2] =-1.0e-5 ;//constant |
| | 302 | fDpi0[3] = 1.50 ; fDpi0[4] = 398. ; fDpi0[5] = 12. ;//mean |
| | 303 | fDpi0[6] =-7.0E-2 ; fDpi0[7] =-524. ; fDpi0[8] = 22. ;//sigma |
| | 304 | fDpi0[9] = 110.; //for E> fDpi0[9] parameters calculated at fDpi0[9] |
| | 305 | |
| | 306 | fDhadron[0] = 6.5 ; fDhadron[1] =-5.3 ; fDhadron[2] = 1.5 ;//constant |
| | 307 | fDhadron[3] = 3.8 ; fDhadron[4] = 0.23 ; fDhadron[5] =-1.2 ;//mean |
| | 308 | fDhadron[6] = 0.88 ; fDhadron[7] = 9.3E-2 ; fDhadron[8] =-0.51 ;//sigma |
| | 309 | fDhadron[9] = 2.; //for E> fDhadron[9] parameters calculated at fDhadron[9] |
| | 310 | |
| | 311 | fDmuon[0] = 0.0631 ; |
| | 312 | fDmuon[1] = 1.4 ; |
| | 313 | fDmuon[2] = 0.0557 ; |
| | 314 | fDFmuon = new TFormula("Shower shape response to muons" , "landau") ; |
| | 315 | fDFmuon->SetParameters( fDmuon[0], fDmuon[1], fDmuon[2]) ; |
| | 316 | |
| | 317 | |
| | 318 | // x(CPV-EMC) distance gaussian parameters |
| | 319 | |
| | 320 | // fXelectron[0] = 8.06e-2 ; fXelectron[1] = 1.00e-2; fXelectron[2] =-5.14e-2;//constant |
| | 321 | // fXelectron[3] = 0.202 ; fXelectron[4] = 8.15e-3; fXelectron[5] = 4.55 ;//mean |
| | 322 | // fXelectron[6] = 0.334 ; fXelectron[7] = 0.186 ; fXelectron[8] = 4.32e-2;//sigma |
| | 323 | |
| | 324 | // //charged hadrons gaus |
| | 325 | // fXcharged[0] = 6.43e-3 ; fXcharged[1] =-4.19e-5; fXcharged[2] = 1.42e-3;//constant |
| | 326 | // fXcharged[3] = 2.75 ; fXcharged[4] =-0.40 ; fXcharged[5] = 1.68 ;//mean |
| | 327 | // fXcharged[6] = 3.135 ; fXcharged[7] =-9.41e-2; fXcharged[8] = 1.31e-2;//sigma |
| | 328 | |
| | 329 | // // z(CPV-EMC) distance gaussian parameters |
| | 330 | |
| | 331 | // fZelectron[0] = 8.22e-2 ; fZelectron[1] = 5.11e-3; fZelectron[2] =-3.05e-2;//constant |
| | 332 | // fZelectron[3] = 3.09e-2 ; fZelectron[4] = 5.87e-2; fZelectron[5] =-9.49e-2;//mean |
| | 333 | // fZelectron[6] = 0.263 ; fZelectron[7] =-9.02e-3; fZelectron[8] = 0.151 ;//sigma |
| | 334 | |
| | 335 | // //charged hadrons gaus |
| | 336 | |
| | 337 | // fZcharged[0] = 1.00e-2 ; fZcharged[1] = 2.82E-4 ; fZcharged[2] = 2.87E-3 ;//constant |
| | 338 | // fZcharged[3] =-4.68e-2 ; fZcharged[4] =-9.21e-3 ; fZcharged[5] = 4.91e-2 ;//mean |
| | 339 | // fZcharged[6] = 1.425 ; fZcharged[7] =-5.90e-2 ; fZcharged[8] = 5.07e-2 ;//sigma |
| | 340 | |
| | 341 | |
| | 342 | fXelectron[0] =-1.6E-2 ; fXelectron[1] = 0.77 ; fXelectron[2] =-0.15 ;//constant |
| | 343 | fXelectron[3] = 0.35 ; fXelectron[4] = 0.25 ; fXelectron[5] = 4.12 ;//mean |
| | 344 | fXelectron[6] = 0.30 ; fXelectron[7] = 0.11 ; fXelectron[8] = 0.16 ;//sigma |
| | 345 | fXelectron[9] = 3.; //for E> fXelectron[9] parameters calculated at fXelectron[9] |
| | 346 | |
| | 347 | //charged hadrons gaus |
| | 348 | fXcharged[0] = 0.14 ; fXcharged[1] =-3.0E-2 ; fXcharged[2] = 0 ;//constant |
| | 349 | fXcharged[3] = 1.4 ; fXcharged[4] =-9.3E-2 ; fXcharged[5] = 1.4 ;//mean |
| | 350 | fXcharged[6] = 5.7 ; fXcharged[7] = 0.27 ; fXcharged[8] =-1.8 ;//sigma |
| | 351 | fXcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] |
| | 352 | |
| | 353 | // z(CPV-EMC) distance gaussian parameters |
| | 354 | |
| | 355 | fZelectron[0] = 0.49 ; fZelectron[1] = 0.53 ; fZelectron[2] =-9.8E-2 ;//constant |
| | 356 | fZelectron[3] = 2.8E-2 ; fZelectron[4] = 5.0E-2 ; fZelectron[5] =-8.2E-2 ;//mean |
| | 357 | fZelectron[6] = 0.25 ; fZelectron[7] =-1.7E-2 ; fZelectron[8] = 0.17 ;//sigma |
| | 358 | fZelectron[9] = 3.; //for E> fZelectron[9] parameters calculated at fZelectron[9] |
| | 359 | |
| | 360 | //charged hadrons gaus |
| | 361 | |
| | 362 | fZcharged[0] = 0.46 ; fZcharged[1] =-0.65 ; fZcharged[2] = 0.52 ;//constant |
| | 363 | fZcharged[3] = 1.1E-2 ; fZcharged[4] = 0. ; fZcharged[5] = 0. ;//mean |
| | 364 | fZcharged[6] = 0.60 ; fZcharged[7] =-8.2E-2 ; fZcharged[8] = 0.45 ;//sigma |
| | 365 | fZcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] |
| | 366 | |
| | 367 | //Threshold to differentiate between charged and neutral |
| | 368 | fChargedNeutralThreshold = 1e-5; |
| | 369 | fTOFEnThreshold = 2; //Maximum energy to use TOF |
| | 370 | fDispEnThreshold = 0.5; //Minimum energy to use shower shape |
| | 371 | fDispMultThreshold = 3; //Minimum multiplicity to use shower shape |
| | 372 | |
| | 373 | //Weight to hadrons recontructed energy |
| | 374 | |
| | 375 | fERecWeightPar[0] = 0.32 ; |
| | 376 | fERecWeightPar[1] = 3.8 ; |
| | 377 | fERecWeightPar[2] = 5.4E-3 ; |
| | 378 | fERecWeightPar[3] = 5.6E-2 ; |
| | 379 | fERecWeight = new TFormula("Weight for hadrons" , "[0]*exp(-x*[1])+[2]*exp(-x*[3])") ; |
| | 380 | fERecWeight ->SetParameters(fERecWeightPar[0],fERecWeightPar[1] ,fERecWeightPar[2] ,fERecWeightPar[3]) ; |
| | 381 | |
| | 382 | |
| | 383 | for (Int_t i =0; i< AliPID::kSPECIESN ; i++) |
| | 384 | fInitPID[i] = 1.; |
| | 385 | |
| | 386 | } |
| | 387 | |
| | 388 | //________________________________________________________________________ |
| | 389 | void AliPHOSPIDv1::Exec(Option_t *option) |
| | 390 | { |
| | 391 | // Steering method to perform particle reconstruction and identification |
| | 392 | // for the event range from fFirstEvent to fLastEvent. |
| | 393 | // This range is optionally set by SetEventRange(). |
| | 394 | // if fLastEvent=-1 (by default), then process events until the end. |
| | 395 | |
| | 396 | if(strstr(option,"tim")) |
| | 397 | gBenchmark->Start("PHOSPID"); |
| | 398 | |
| | 399 | if(strstr(option,"print")) { |
| | 400 | Print() ; |
| | 401 | return ; |
| | 402 | } |
| | 403 | |
| | 404 | |
| | 405 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
| | 406 | |
| | 407 | if (fLastEvent == -1) |
| | 408 | fLastEvent = gime->MaxEvent() - 1 ; |
| | 409 | else |
| | 410 | fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent()); |
| | 411 | Int_t nEvents = fLastEvent - fFirstEvent + 1; |
| | 412 | |
| | 413 | Int_t ievent ; |
| | 414 | for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { |
| | 415 | gime->Event(ievent,"TR") ; |
| | 416 | if(gime->TrackSegments() && //Skip events, where no track segments made |
| | 417 | gime->TrackSegments()->GetEntriesFast()) { |
| | 418 | |
| | 419 | MakeRecParticles() ; |
| | 420 | |
| | 421 | if(fBayesian) |
| | 422 | MakePID() ; |
| | 423 | |
| | 424 | WriteRecParticles(); |
| | 425 | if(strstr(option,"deb")) |
| | 426 | PrintRecParticles(option) ; |
| | 427 | //increment the total number of rec particles per run |
| | 428 | fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; |
| | 429 | } |
| | 430 | } |
| | 431 | if(strstr(option,"deb")) |
| | 432 | PrintRecParticles(option); |
| | 433 | if(strstr(option,"tim")){ |
| | 434 | gBenchmark->Stop("PHOSPID"); |
| | 435 | AliInfo(Form("took %f seconds for PID %f seconds per event", |
| | 436 | gBenchmark->GetCpuTime("PHOSPID"), |
| | 437 | gBenchmark->GetCpuTime("PHOSPID")/nEvents)) ; |
| | 438 | } |
| | 439 | if(fWrite) |
| | 440 | Unload(); |
| | 441 | } |
| | 442 | |
| | 443 | //________________________________________________________________________ |
| | 444 | Double_t AliPHOSPIDv1::GausF(Double_t x, Double_t y, Double_t * par) |
| | 445 | { |
| | 446 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
| | 447 | //this method returns a density probability of this parameter, given by a gaussian |
| | 448 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b |
| | 449 | //Float_t xorg = x; |
| | 450 | if (x > par[9]) x = par[9]; |
| | 451 | |
| | 452 | //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; |
| | 453 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
| | 454 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
| | 455 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; |
| | 456 | |
| | 457 | // if(xorg > 30) |
| | 458 | // cout<<"En_in = "<<xorg<<"; En_out = "<<x<<"; cnt = "<<cnt |
| | 459 | // <<"; mean = "<<mean<<"; sigma = "<<sigma<<endl; |
| | 460 | |
| | 461 | // Double_t arg = - (y-mean) * (y-mean) / (2*sigma*sigma) ; |
| | 462 | // return cnt * TMath::Exp(arg) ; |
| | 463 | if(TMath::Abs(sigma) > 1.e-10){ |
| | 464 | return cnt*TMath::Gaus(y,mean,sigma); |
| | 465 | } |
| | 466 | else |
| | 467 | return 0.; |
| | 468 | |
| | 469 | } |
| | 470 | //________________________________________________________________________ |
| | 471 | Double_t AliPHOSPIDv1::GausPol2(Double_t x, Double_t y, Double_t * par) |
| | 472 | { |
| | 473 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
| | 474 | //this method returns a density probability of this parameter, given by a gaussian |
| | 475 | //function whose parameters depend with the energy like second order polinomial |
| | 476 | |
| | 477 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
| | 478 | Double_t mean = par[3] + par[4] * x + par[5] * x * x ; |
| | 479 | Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; |
| | 480 | |
| | 481 | if(TMath::Abs(sigma) > 1.e-10){ |
| | 482 | return cnt*TMath::Gaus(y,mean,sigma); |
| | 483 | } |
| | 484 | else |
| | 485 | return 0.; |
| | 486 | |
| | 487 | |
| | 488 | |
| | 489 | } |
| | 490 | |
| | 491 | //____________________________________________________________________________ |
| | 492 | const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const |
| | 493 | { |
| | 494 | //Get file name that contains the PCA for a particle ("photon or pi0") |
| | 495 | particle.ToLower(); |
| | 496 | TString name; |
| | 497 | if (particle=="photon") |
| | 498 | name = fFileNamePrincipalPhoton ; |
| | 499 | else if (particle=="pi0" ) |
| | 500 | name = fFileNamePrincipalPi0 ; |
| | 501 | else |
| | 502 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
| | 503 | particle.Data())); |
| | 504 | return name; |
| | 505 | } |
| | 506 | |
| | 507 | //____________________________________________________________________________ |
| | 508 | Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const |
| | 509 | { |
| | 510 | // Get the i-th parameter "Calibration" |
| | 511 | Float_t param = 0.; |
| | 512 | if (i>2 || i<0) { |
| | 513 | AliError(Form("Invalid parameter number: %d",i)); |
| | 514 | } else |
| | 515 | param = (*fParameters)(0,i); |
| | 516 | return param; |
| | 517 | } |
| | 518 | |
| | 519 | //____________________________________________________________________________ |
| | 520 | Float_t AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const |
| | 521 | { |
| | 522 | // It calibrates Energy depending on the recpoint energy. |
| | 523 | // The energy of the reconstructed cluster is corrected with |
| | 524 | // the formula A + B* E + C* E^2, whose parameters where obtained |
| | 525 | // through the study of the reconstructed energy distribution of |
| | 526 | // monoenergetic photons. |
| | 527 | |
| | 528 | Float_t p[]={0.,0.,0.}; |
| | 529 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); |
| | 530 | Float_t enerec = p[0] + p[1]*e + p[2]*e*e; |
| | 531 | return enerec ; |
| | 532 | |
| | 533 | } |
| | 534 | |
| | 535 | //____________________________________________________________________________ |
| | 536 | Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const |
| | 537 | { |
| | 538 | // Get the i-th parameter "CPV-EMC distance" for the specified axis |
| | 539 | Float_t param = 0.; |
| | 540 | if(i>2 || i<0) { |
| | 541 | AliError(Form("Invalid parameter number: %d",i)); |
| | 542 | } else { |
| | 543 | axis.ToLower(); |
| | 544 | if (axis == "x") |
| | 545 | param = (*fParameters)(1,i); |
| | 546 | else if (axis == "z") |
| | 547 | param = (*fParameters)(2,i); |
| | 548 | else { |
| | 549 | AliError(Form("Invalid axis name: %s",axis.Data())); |
| | 550 | } |
| | 551 | } |
| | 552 | return param; |
| | 553 | } |
| | 554 | |
| | 555 | //____________________________________________________________________________ |
| | 556 | Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const |
| | 557 | { |
| | 558 | // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and |
| | 559 | // Purity-Efficiency point |
| | 560 | |
| | 561 | axis.ToLower(); |
| | 562 | Float_t p[]={0.,0.,0.}; |
| | 563 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis); |
| | 564 | Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e); |
| | 565 | return sig; |
| | 566 | } |
| | 567 | |
| | 568 | //____________________________________________________________________________ |
| | 569 | Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const |
| | 570 | { |
| | 571 | // Calculates the parameter param of the ellipse |
| | 572 | |
| | 573 | particle.ToLower(); |
| | 574 | param. ToLower(); |
| | 575 | Float_t p[4]={0.,0.,0.,0.}; |
| | 576 | Float_t value = 0.0; |
| | 577 | for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i); |
| | 578 | if (particle == "photon") { |
| | 579 | if (param.Contains("a")) e = TMath::Min((Double_t)e,70.); |
| | 580 | else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.); |
| | 581 | else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1); |
| | 582 | } |
| | 583 | |
| | 584 | if (particle == "photon") |
| | 585 | value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3]; |
| | 586 | else if (particle == "pi0") |
| | 587 | value = p[0] + p[1]*e + p[2]*e*e; |
| | 588 | |
| | 589 | return value; |
| | 590 | } |
| | 591 | |
| | 592 | //_____________________________________________________________________________ |
| | 593 | Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const |
| | 594 | { |
| | 595 | // Get the parameter "i" to calculate the boundary on the moment M2x |
| | 596 | // for photons at high p_T |
| | 597 | Float_t param = 0; |
| | 598 | if (i>3 || i<0) { |
| | 599 | AliError(Form("Wrong parameter number: %d\n",i)); |
| | 600 | } else |
| | 601 | param = (*fParameters)(14,i) ; |
| | 602 | return param; |
| | 603 | } |
| | 604 | |
| | 605 | //____________________________________________________________________________ |
| | 606 | Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const |
| | 607 | { |
| | 608 | // Get the parameter "i" to calculate the boundary on the moment M2x |
| | 609 | // for pi0 at high p_T |
| | 610 | Float_t param = 0; |
| | 611 | if (i>2 || i<0) { |
| | 612 | AliError(Form("Wrong parameter number: %d\n",i)); |
| | 613 | } else |
| | 614 | param = (*fParameters)(15,i) ; |
| | 615 | return param; |
| | 616 | } |
| | 617 | |
| | 618 | //____________________________________________________________________________ |
| | 619 | Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const |
| | 620 | { |
| | 621 | // Get TimeGate parameter depending on Purity-Efficiency i: |
| | 622 | // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity |
| | 623 | Float_t param = 0.; |
| | 624 | if(i>2 || i<0) { |
| | 625 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); |
| | 626 | } else |
| | 627 | param = (*fParameters)(3,i) ; |
| | 628 | return param; |
| | 629 | } |
| | 630 | |
| | 631 | //_____________________________________________________________________________ |
| | 632 | Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const |
| | 633 | { |
| | 634 | // Get the parameter "i" that is needed to calculate the ellipse |
| | 635 | // parameter "param" for the particle "particle" ("photon" or "pi0") |
| | 636 | |
| | 637 | particle.ToLower(); |
| | 638 | param. ToLower(); |
| | 639 | Int_t offset = -1; |
| | 640 | if (particle == "photon") |
| | 641 | offset=0; |
| | 642 | else if (particle == "pi0") |
| | 643 | offset=5; |
| | 644 | else |
| | 645 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
| | 646 | particle.Data())); |
| | 647 | |
| | 648 | Int_t p= -1; |
| | 649 | Float_t par = 0; |
| | 650 | |
| | 651 | if (param.Contains("a")) p=4+offset; |
| | 652 | else if(param.Contains("b")) p=5+offset; |
| | 653 | else if(param.Contains("c")) p=6+offset; |
| | 654 | else if(param.Contains("x0"))p=7+offset; |
| | 655 | else if(param.Contains("y0"))p=8+offset; |
| | 656 | |
| | 657 | if (i>4 || i<0) { |
| | 658 | AliError(Form("No parameter with index %d", i)) ; |
| | 659 | } else if (p==-1) { |
| | 660 | AliError(Form("No parameter with name %s", param.Data() )) ; |
| | 661 | } else |
| | 662 | par = (*fParameters)(p,i) ; |
| | 663 | |
| | 664 | return par; |
| | 665 | } |
| | 666 | |
| | 667 | |
| | 668 | //____________________________________________________________________________ |
| | 669 | Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const |
| | 670 | { |
| | 671 | // Calculates the distance between the EMC RecPoint and the PPSD RecPoint |
| | 672 | |
| | 673 | const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
| | 674 | TVector3 vecEmc ; |
| | 675 | TVector3 vecCpv ; |
| | 676 | if(cpv){ |
| | 677 | emc->GetLocalPosition(vecEmc) ; |
| | 678 | cpv->GetLocalPosition(vecCpv) ; |
| | 679 | |
| | 680 | if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ |
| | 681 | // Correct to difference in CPV and EMC position due to different distance to center. |
| | 682 | // we assume, that particle moves from center |
| | 683 | Float_t dCPV = geom->GetIPtoOuterCoverDistance(); |
| | 684 | Float_t dEMC = geom->GetIPtoCrystalSurface() ; |
| | 685 | dEMC = dEMC / dCPV ; |
| | 686 | vecCpv = dEMC * vecCpv - vecEmc ; |
| | 687 | if (axis == "X") return vecCpv.X(); |
| | 688 | if (axis == "Y") return vecCpv.Y(); |
| | 689 | if (axis == "Z") return vecCpv.Z(); |
| | 690 | if (axis == "R") return vecCpv.Mag(); |
| | 691 | } |
| | 692 | return 100000000 ; |
| | 693 | } |
| | 694 | return 100000000 ; |
| | 695 | } |
| | 696 | //____________________________________________________________________________ |
| | 697 | Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Int_t effPur, Float_t e) const |
| | 698 | { |
| | 699 | //Calculates the pid bit for the CPV selection per each purity. |
| | 700 | if(effPur>2 || effPur<0) |
| | 701 | AliError(Form("Invalid Efficiency-Purity choice %d",effPur)); |
| | 702 | |
| | 703 | Float_t sigX = GetCpv2EmcDistanceCut("X",e); |
| | 704 | Float_t sigZ = GetCpv2EmcDistanceCut("Z",e); |
| | 705 | |
| | 706 | Float_t deltaX = TMath::Abs(GetDistance(emc, cpv, "X")); |
| | 707 | Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv, "Z")); |
| | 708 | //Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ; |
| | 709 | |
| | 710 | //if(deltaX>sigX*(effPur+1)) |
| | 711 | //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1))) |
| | 712 | if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1))) |
| | 713 | return 1;//Neutral |
| | 714 | else |
| | 715 | return 0;//Charged |
| | 716 | } |
| | 717 | |
| | 718 | //____________________________________________________________________________ |
| | 719 | Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const |
| | 720 | { |
| | 721 | //Is the particle inside de PCA ellipse? |
| | 722 | |
| | 723 | particle.ToLower(); |
| | 724 | Int_t prinbit = 0 ; |
| | 725 | Float_t a = GetEllipseParameter(particle,"a" , e); |
| | 726 | Float_t b = GetEllipseParameter(particle,"b" , e); |
| | 727 | Float_t c = GetEllipseParameter(particle,"c" , e); |
| | 728 | Float_t x0 = GetEllipseParameter(particle,"x0", e); |
| | 729 | Float_t y0 = GetEllipseParameter(particle,"y0", e); |
| | 730 | |
| | 731 | Float_t r = TMath::Power((p[0] - x0)/a,2) + |
| | 732 | TMath::Power((p[1] - y0)/b,2) + |
| | 733 | c*(p[0] - x0)*(p[1] - y0)/(a*b) ; |
| | 734 | //3 different ellipses defined |
| | 735 | if((effPur==2) && (r<1./2.)) prinbit= 1; |
| | 736 | if((effPur==1) && (r<2. )) prinbit= 1; |
| | 737 | if((effPur==0) && (r<9./2.)) prinbit= 1; |
| | 738 | |
| | 739 | if(r<0) |
| | 740 | AliError("Negative square?") ; |
| | 741 | |
| | 742 | return prinbit; |
| | 743 | |
| | 744 | } |
| | 745 | //____________________________________________________________________________ |
| | 746 | Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const |
| | 747 | { |
| | 748 | // Set bit for identified hard photons (E > 30 GeV) |
| | 749 | // if the second moment M2x is below the boundary |
| | 750 | |
| | 751 | Float_t e = emc->GetEnergy(); |
| | 752 | if (e < 30.0) return 0; |
| | 753 | Float_t m2x = emc->GetM2x(); |
| | 754 | Float_t m2xBoundary = GetParameterPhotonBoundary(0) * |
| | 755 | TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ |
| | 756 | TMath::Power(GetParameterPhotonBoundary(2),2)) + |
| | 757 | GetParameterPhotonBoundary(3); |
| | 758 | AliDebug(1, Form("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f", |
| | 759 | e,m2x,m2xBoundary)); |
| | 760 | if (m2x < m2xBoundary) |
| | 761 | return 1;// A hard photon |
| | 762 | else |
| | 763 | return 0;// Not a hard photon |
| | 764 | } |
| | 765 | |
| | 766 | //____________________________________________________________________________ |
| | 767 | Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const |
| | 768 | { |
| | 769 | // Set bit for identified hard pi0 (E > 30 GeV) |
| | 770 | // if the second moment M2x is above the boundary |
| | 771 | |
| | 772 | Float_t e = emc->GetEnergy(); |
| | 773 | if (e < 30.0) return 0; |
| | 774 | Float_t m2x = emc->GetM2x(); |
| | 775 | Float_t m2xBoundary = GetParameterPi0Boundary(0) + |
| | 776 | e * GetParameterPi0Boundary(1); |
| | 777 | AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary)); |
| | 778 | if (m2x > m2xBoundary) |
| | 779 | return 1;// A hard pi0 |
| | 780 | else |
| | 781 | return 0;// Not a hard pi0 |
| | 782 | } |
| | 783 | |
| | 784 | //____________________________________________________________________________ |
| | 785 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const |
| | 786 | { |
| | 787 | // Calculates the momentum direction: |
| | 788 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint |
| | 789 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints |
| | 790 | // However because of the poor position resolution of PPSD the direction is always taken as if we were |
| | 791 | // in case 1. |
| | 792 | |
| | 793 | TVector3 dir(0,0,0) ; |
| | 794 | TMatrixF dummy ; |
| | 795 | |
| | 796 | emc->GetGlobalPosition(dir, dummy) ; |
| | 797 | |
| | 798 | //account correction to the position of IP |
| | 799 | Float_t xo,yo,zo ; //Coordinates of the origin |
| | 800 | if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){ |
| | 801 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; |
| | 802 | } |
| | 803 | else{ |
| | 804 | xo=yo=zo=0.; |
| | 805 | } |
| | 806 | TVector3 origin(xo,yo,zo); |
| | 807 | dir = dir - origin ; |
| | 808 | dir.SetMag(1.) ; |
| | 809 | |
| | 810 | return dir ; |
| | 811 | } |
| | 812 | |
| | 813 | //________________________________________________________________________ |
| | 814 | Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par) |
| | 815 | { |
| | 816 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
| | 817 | //this method returns a density probability of this parameter, given by a landau |
| | 818 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b |
| | 819 | |
| | 820 | if (x > par[9]) x = par[9]; |
| | 821 | |
| | 822 | //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; |
| | 823 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
| | 824 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
| | 825 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; |
| | 826 | |
| | 827 | if(TMath::Abs(sigma) > 1.e-10){ |
| | 828 | return cnt*TMath::Landau(y,mean,sigma); |
| | 829 | } |
| | 830 | else |
| | 831 | return 0.; |
| | 832 | |
| | 833 | } |
| | 834 | //________________________________________________________________________ |
| | 835 | Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par) |
| | 836 | { |
| | 837 | |
| | 838 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
| | 839 | //this method returns a density probability of this parameter, given by a landau |
| | 840 | //function whose parameters depend with the energy like second order polinomial |
| | 841 | |
| | 842 | Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ; |
| | 843 | Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ; |
| | 844 | Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ; |
| | 845 | |
| | 846 | if(TMath::Abs(sigma) > 1.e-10){ |
| | 847 | return cnt*TMath::Landau(y,mean,sigma); |
| | 848 | } |
| | 849 | else |
| | 850 | return 0.; |
| | 851 | |
| | 852 | |
| | 853 | } |
| | 854 | // //________________________________________________________________________ |
| | 855 | // Double_t AliPHOSPIDv1::ChargedHadronDistProb(Double_t x, Double_t y, Double_t * parg, Double_t * parl) |
| | 856 | // { |
| | 857 | // Double_t cnt = 0.0 ; |
| | 858 | // Double_t mean = 0.0 ; |
| | 859 | // Double_t sigma = 0.0 ; |
| | 860 | // Double_t arg = 0.0 ; |
| | 861 | // if (y < parl[4] / (x*x) + parl[5] / x + parl[3]){ |
| | 862 | // cnt = parg[1] / (x*x) + parg[2] / x + parg[0] ; |
| | 863 | // mean = parg[4] / (x*x) + parg[5] / x + parg[3] ; |
| | 864 | // sigma = parg[7] / (x*x) + parg[8] / x + parg[6] ; |
| | 865 | // TF1 * f = new TF1("gaus","gaus",0.,100.); |
| | 866 | // f->SetParameters(cnt,mean,sigma); |
| | 867 | // arg = f->Eval(y) ; |
| | 868 | // } |
| | 869 | // else{ |
| | 870 | // cnt = parl[1] / (x*x) + parl[2] / x + parl[0] ; |
| | 871 | // mean = parl[4] / (x*x) + parl[5] / x + parl[3] ; |
| | 872 | // sigma = parl[7] / (x*x) + parl[8] / x + parl[6] ; |
| | 873 | // TF1 * f = new TF1("landau","landau",0.,100.); |
| | 874 | // f->SetParameters(cnt,mean,sigma); |
| | 875 | // arg = f->Eval(y) ; |
| | 876 | // } |
| | 877 | // // Double_t mean = par[3] + par[4] * x + par[5] * x * x ; |
| | 878 | // // Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; |
| | 879 | |
| | 880 | // //Double_t arg = -(y-mean)*(y-mean)/(2*sigma*sigma) ; |
| | 881 | // //return cnt * TMath::Exp(arg) ; |
| | 882 | |
| | 883 | // return arg; |
| | 884 | |
| | 885 | // } |
| | 886 | //____________________________________________________________________________ |
| | 887 | void AliPHOSPIDv1::MakePID() |
| | 888 | { |
| | 889 | // construct the PID weight from a Bayesian Method |
| | 890 | |
| | 891 | const Int_t kSPECIES = AliPID::kSPECIESN ; |
| | 892 | |
| | 893 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
| | 894 | |
| | 895 | Int_t nparticles = gime->RecParticles()->GetEntriesFast() ; |
| | 896 | |
| | 897 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
| | 898 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
| | 899 | TClonesArray * trackSegments = gime->TrackSegments() ; |
| | 900 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { |
| | 901 | AliFatal("RecPoints or TrackSegments not found !") ; |
| | 902 | } |
| | 903 | TIter next(trackSegments) ; |
| | 904 | AliPHOSTrackSegment * ts ; |
| | 905 | Int_t index = 0 ; |
| | 906 | |
| | 907 | Double_t * stof[kSPECIES] ; |
| | 908 | Double_t * sdp [kSPECIES] ; |
| | 909 | Double_t * scpv[kSPECIES] ; |
| | 910 | Double_t * sw [kSPECIES] ; |
| | 911 | //Info("MakePID","Begin MakePID"); |
| | 912 | |
| | 913 | for (Int_t i =0; i< kSPECIES; i++){ |
| | 914 | stof[i] = new Double_t[nparticles] ; |
| | 915 | sdp [i] = new Double_t[nparticles] ; |
| | 916 | scpv[i] = new Double_t[nparticles] ; |
| | 917 | sw [i] = new Double_t[nparticles] ; |
| | 918 | } |
| | 919 | |
| | 920 | |
| | 921 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
| | 922 | |
| | 923 | //cout<<">>>>>> Bayesian Index "<<index<<endl; |
| | 924 | |
| | 925 | AliPHOSEmcRecPoint * emc = 0 ; |
| | 926 | if(ts->GetEmcIndex()>=0) |
| | 927 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; |
| | 928 | |
| | 929 | AliPHOSCpvRecPoint * cpv = 0 ; |
| | 930 | if(ts->GetCpvIndex()>=0) |
| | 931 | cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
| | 932 | |
| | 933 | // Int_t track = 0 ; |
| | 934 | // track = ts->GetTrackIndex() ; //TPC tracks ? |
| | 935 | |
| | 936 | if (!emc) { |
| | 937 | AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; |
| | 938 | } |
| | 939 | |
| | 940 | |
| | 941 | // ############Tof############################# |
| | 942 | |
| | 943 | // Info("MakePID", "TOF"); |
| | 944 | Float_t en = emc->GetEnergy(); |
| | 945 | Double_t time = emc->GetTime() ; |
| | 946 | // cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl; |
| | 947 | |
| | 948 | // now get the signals probability |
| | 949 | // s(pid) in the Bayesian formulation |
| | 950 | |
| | 951 | stof[AliPID::kPhoton][index] = 1.; |
| | 952 | stof[AliPID::kElectron][index] = 1.; |
| | 953 | stof[AliPID::kEleCon][index] = 1.; |
| | 954 | //We assing the same prob to charged hadrons, sum is 1 |
| | 955 | stof[AliPID::kPion][index] = 1./3.; |
| | 956 | stof[AliPID::kKaon][index] = 1./3.; |
| | 957 | stof[AliPID::kProton][index] = 1./3.; |
| | 958 | //We assing the same prob to neutral hadrons, sum is 1 |
| | 959 | stof[AliPID::kNeutron][index] = 1./2.; |
| | 960 | stof[AliPID::kKaon0][index] = 1./2.; |
| | 961 | stof[AliPID::kMuon][index] = 1.; |
| | 962 | |
| | 963 | if(en < fTOFEnThreshold) { |
| | 964 | |
| | 965 | Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution |
| | 966 | Double_t pTofKaon = 0; |
| | 967 | |
| | 968 | if(time < fTkaonl[1]) |
| | 969 | pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution |
| | 970 | else |
| | 971 | pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution |
| | 972 | |
| | 973 | Double_t pTofNucleon = 0; |
| | 974 | |
| | 975 | if(time < fThhadronl[1]) |
| | 976 | pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution |
| | 977 | else |
| | 978 | pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution |
| | 979 | //We assing the same prob to neutral hadrons, sum is the average prob |
| | 980 | Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ; |
| | 981 | //We assing the same prob to charged hadrons, sum is the average prob |
| | 982 | Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ; |
| | 983 | |
| | 984 | stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ; |
| | 985 | //gaus distribution |
| | 986 | stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ; |
| | 987 | //a conversion electron has the photon ToF |
| | 988 | stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ; |
| | 989 | |
| | 990 | stof[AliPID::kElectron][index] = pTofPion ; |
| | 991 | |
| | 992 | stof[AliPID::kPion][index] = pTofChHadron ; |
| | 993 | stof[AliPID::kKaon][index] = pTofChHadron ; |
| | 994 | stof[AliPID::kProton][index] = pTofChHadron ; |
| | 995 | |
| | 996 | stof[AliPID::kKaon0][index] = pTofNeHadron ; |
| | 997 | stof[AliPID::kNeutron][index] = pTofNeHadron ; |
| | 998 | } |
| | 999 | |
| | 1000 | // Info("MakePID", "Dispersion"); |
| | 1001 | |
| | 1002 | // ###########Shower shape: Dispersion#################### |
| | 1003 | Float_t dispersion = emc->GetDispersion(); |
| | 1004 | //dispersion is not well defined if the cluster is only in few crystals |
| | 1005 | |
| | 1006 | sdp[AliPID::kPhoton][index] = 1. ; |
| | 1007 | sdp[AliPID::kElectron][index] = 1. ; |
| | 1008 | sdp[AliPID::kPion][index] = 1. ; |
| | 1009 | sdp[AliPID::kKaon][index] = 1. ; |
| | 1010 | sdp[AliPID::kProton][index] = 1. ; |
| | 1011 | sdp[AliPID::kNeutron][index] = 1. ; |
| | 1012 | sdp[AliPID::kEleCon][index] = 1. ; |
| | 1013 | sdp[AliPID::kKaon0][index] = 1. ; |
| | 1014 | sdp[AliPID::kMuon][index] = 1. ; |
| | 1015 | |
| | 1016 | if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){ |
| | 1017 | sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ; |
| | 1018 | sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ; |
| | 1019 | sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ; |
| | 1020 | sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ; |
| | 1021 | sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ; |
| | 1022 | sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ; |
| | 1023 | sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index]; |
| | 1024 | sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ; |
| | 1025 | sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ; |
| | 1026 | //landau distribution |
| | 1027 | } |
| | 1028 | |
| | 1029 | // Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion); |
| | 1030 | // Info("MakePID","ss: photon %f, hadron %f ", sdp[AliPID::kPhoton][index], sdp[AliPID::kPion][index]); |
| | 1031 | // cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ; |
| | 1032 | // cout<<"<<<<<ss: photon "<<sdp[AliPID::kPhoton][index]<<", hadron "<<sdp[AliPID::kPion][index]<<endl; |
| | 1033 | |
| | 1034 | //########## CPV-EMC Distance####################### |
| | 1035 | // Info("MakePID", "Distance"); |
| | 1036 | |
| | 1037 | Float_t x = TMath::Abs(GetDistance(emc, cpv, "X")) ; |
| | 1038 | Float_t z = GetDistance(emc, cpv, "Z") ; |
| | 1039 | |
| | 1040 | Double_t pcpv = 0 ; |
| | 1041 | Double_t pcpvneutral = 0. ; |
| | 1042 | |
| | 1043 | Double_t elprobx = GausF(en , x, fXelectron) ; |
| | 1044 | Double_t elprobz = GausF(en , z, fZelectron) ; |
| | 1045 | Double_t chprobx = GausF(en , x, fXcharged) ; |
| | 1046 | Double_t chprobz = GausF(en , z, fZcharged) ; |
| | 1047 | Double_t pcpvelectron = elprobx * elprobz; |
| | 1048 | Double_t pcpvcharged = chprobx * chprobz; |
| | 1049 | |
| | 1050 | // cout<<">>>>energy "<<en<<endl; |
| | 1051 | // cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl; |
| | 1052 | // cout<<">>>>hadron : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl; |
| | 1053 | // cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; |
| | 1054 | |
| | 1055 | // Is neutral or charged? |
| | 1056 | if(pcpvelectron >= pcpvcharged) |
| | 1057 | pcpv = pcpvelectron ; |
| | 1058 | else |
| | 1059 | pcpv = pcpvcharged ; |
| | 1060 | |
| | 1061 | if(pcpv < fChargedNeutralThreshold) |
| | 1062 | { |
| | 1063 | pcpvneutral = 1. ; |
| | 1064 | pcpvcharged = 0. ; |
| | 1065 | pcpvelectron = 0. ; |
| | 1066 | } |
| | 1067 | // else |
| | 1068 | // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl; |
| | 1069 | |
| | 1070 | scpv[AliPID::kPion][index] = pcpvcharged ; |
| | 1071 | scpv[AliPID::kKaon][index] = pcpvcharged ; |
| | 1072 | scpv[AliPID::kProton][index] = pcpvcharged ; |
| | 1073 | |
| | 1074 | scpv[AliPID::kMuon][index] = pcpvelectron ; |
| | 1075 | scpv[AliPID::kElectron][index] = pcpvelectron ; |
| | 1076 | scpv[AliPID::kEleCon][index] = pcpvelectron ; |
| | 1077 | |
| | 1078 | scpv[AliPID::kPhoton][index] = pcpvneutral ; |
| | 1079 | scpv[AliPID::kNeutron][index] = pcpvneutral ; |
| | 1080 | scpv[AliPID::kKaon0][index] = pcpvneutral ; |
| | 1081 | |
| | 1082 | |
| | 1083 | // Info("MakePID", "CPV passed"); |
| | 1084 | |
| | 1085 | //############## Pi0 ############################# |
| | 1086 | stof[AliPID::kPi0][index] = 0. ; |
| | 1087 | scpv[AliPID::kPi0][index] = 0. ; |
| | 1088 | sdp [AliPID::kPi0][index] = 0. ; |
| | 1089 | |
| | 1090 | if(en > 30.){ |
| | 1091 | // pi0 are detected via decay photon |
| | 1092 | stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index]; |
| | 1093 | scpv[AliPID::kPi0][index] = pcpvneutral ; |
| | 1094 | if(emc->GetMultiplicity() > fDispMultThreshold) |
| | 1095 | sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ; |
| | 1096 | //sdp [AliPID::kPi0][index] = GausPol2(en , dispersion, fDpi0) ; |
| | 1097 | // cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = " |
| | 1098 | // <<emc->GetMultiplicity()<<endl; |
| | 1099 | // cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = " |
| | 1100 | // <<sdp [AliPID::kPi0][index]<<endl; |
| | 1101 | } |
| | 1102 | |
| | 1103 | |
| | 1104 | |
| | 1105 | |
| | 1106 | //############## muon ############################# |
| | 1107 | |
| | 1108 | if(en > 0.5){ |
| | 1109 | //Muons deposit few energy |
| | 1110 | scpv[AliPID::kMuon][index] = 0 ; |
| | 1111 | stof[AliPID::kMuon][index] = 0 ; |
| | 1112 | sdp [AliPID::kMuon][index] = 0 ; |
| | 1113 | } |
| | 1114 | |
| | 1115 | //Weight to apply to hadrons due to energy reconstruction |
| | 1116 | |
| | 1117 | Float_t weight = fERecWeight ->Eval(en) ; |
| | 1118 | |
| | 1119 | sw[AliPID::kPhoton][index] = 1. ; |
| | 1120 | sw[AliPID::kElectron][index] = 1. ; |
| | 1121 | sw[AliPID::kPion][index] = weight ; |
| | 1122 | sw[AliPID::kKaon][index] = weight ; |
| | 1123 | sw[AliPID::kProton][index] = weight ; |
| | 1124 | sw[AliPID::kNeutron][index] = weight ; |
| | 1125 | sw[AliPID::kEleCon][index] = 1. ; |
| | 1126 | sw[AliPID::kKaon0][index] = weight ; |
| | 1127 | sw[AliPID::kMuon][index] = weight ; |
| | 1128 | sw[AliPID::kPi0][index] = 1. ; |
| | 1129 | |
| | 1130 | // if(en > 0.5){ |
| | 1131 | // cout<<"######################################################"<<endl; |
| | 1132 | // //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ; |
| | 1133 | // cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ; |
| | 1134 | // cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl; |
| | 1135 | // cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl; |
| | 1136 | // cout<<">>>>hadron : xprob "<<chprobx<<" zprob "<<chprobz<<endl; |
| | 1137 | // cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; |
| | 1138 | |
| | 1139 | // cout<<"Photon , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index] |
| | 1140 | // <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl; |
| | 1141 | // cout<<"EleCon , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index] |
| | 1142 | // <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl; |
| | 1143 | // cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index] |
| | 1144 | // <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl; |
| | 1145 | // cout<<"Muon , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index] |
| | 1146 | // <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl; |
| | 1147 | // cout<<"Pi0 , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index] |
| | 1148 | // <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl; |
| | 1149 | // cout<<"Pion , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index] |
| | 1150 | // <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl; |
| | 1151 | // cout<<"Kaon0 , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index] |
| | 1152 | // <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl; |
| | 1153 | // cout<<"Kaon , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index] |
| | 1154 | // <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl; |
| | 1155 | // cout<<"Neutron , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index] |
| | 1156 | // <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl; |
| | 1157 | // cout<<"Proton , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index] |
| | 1158 | // <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl; |
| | 1159 | // cout<<"######################################################"<<endl; |
| | 1160 | // } |
| | 1161 | index++; |
| | 1162 | } |
| | 1163 | |
| | 1164 | //for (index = 0 ; index < kSPECIES ; index++) |
| | 1165 | // pid[index] /= nparticles ; |
| | 1166 | |
| | 1167 | |
| | 1168 | // Info("MakePID", "Total Probability calculation"); |
| | 1169 | |
| | 1170 | for(index = 0 ; index < nparticles ; index ++) { |
| | 1171 | |
| | 1172 | AliPHOSRecParticle * recpar = gime->RecParticle(index) ; |
| | 1173 | |
| | 1174 | //Conversion electron? |
| | 1175 | |
| | 1176 | if(recpar->IsEleCon()){ |
| | 1177 | fInitPID[AliPID::kEleCon] = 1. ; |
| | 1178 | fInitPID[AliPID::kPhoton] = 0. ; |
| | 1179 | fInitPID[AliPID::kElectron] = 0. ; |
| | 1180 | } |
| | 1181 | else{ |
| | 1182 | fInitPID[AliPID::kEleCon] = 0. ; |
| | 1183 | fInitPID[AliPID::kPhoton] = 1. ; |
| | 1184 | fInitPID[AliPID::kElectron] = 1. ; |
| | 1185 | } |
| | 1186 | // fInitPID[AliPID::kEleCon] = 0. ; |
| | 1187 | |
| | 1188 | |
| | 1189 | // calculates the Bayesian weight |
| | 1190 | |
| | 1191 | Int_t jndex ; |
| | 1192 | Double_t wn = 0.0 ; |
| | 1193 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) |
| | 1194 | wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * |
| | 1195 | sw[jndex][index] * fInitPID[jndex] ; |
| | 1196 | |
| | 1197 | // cout<<"*************wn "<<wn<<endl; |
| | 1198 | if (TMath::Abs(wn)>0) |
| | 1199 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) { |
| | 1200 | //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn" |
| | 1201 | //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex] << endl; |
| | 1202 | //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl; |
| | 1203 | // if(jndex == AliPID::kPi0 || jndex == AliPID::kPhoton){ |
| | 1204 | // cout<<"Particle "<<jndex<<" final prob * wn " |
| | 1205 | // <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * |
| | 1206 | // fInitPID[jndex] <<" wn "<< wn<<endl; |
| | 1207 | // cout<<"pid "<< fInitPID[jndex]<<", tof "<<stof[jndex][index] |
| | 1208 | // <<", cpv "<<scpv[jndex][index]<<" ss "<<sdp[jndex][index]<<endl; |
| | 1209 | // } |
| | 1210 | recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] * |
| | 1211 | sw[jndex][index] * scpv[jndex][index] * |
| | 1212 | fInitPID[jndex] / wn) ; |
| | 1213 | } |
| | 1214 | } |
| | 1215 | // Info("MakePID", "Delete"); |
| | 1216 | |
| | 1217 | for (Int_t i =0; i< kSPECIES; i++){ |
| | 1218 | delete [] stof[i]; |
| | 1219 | delete [] sdp [i]; |
| | 1220 | delete [] scpv[i]; |
| | 1221 | delete [] sw [i]; |
| | 1222 | } |
| | 1223 | // Info("MakePID","End MakePID"); |
| | 1224 | } |
| | 1225 | |
| | 1226 | //____________________________________________________________________________ |
| | 1227 | void AliPHOSPIDv1::MakeRecParticles() |
| | 1228 | { |
| | 1229 | // Makes a RecParticle out of a TrackSegment |
| | 1230 | |
| | 1231 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
| | 1232 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
| | 1233 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
| | 1234 | TClonesArray * trackSegments = gime->TrackSegments() ; |
| | 1235 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { |
| | 1236 | AliFatal("RecPoints or TrackSegments not found !") ; |
| | 1237 | } |
| | 1238 | TClonesArray * recParticles = gime->RecParticles() ; |
| | 1239 | recParticles->Clear(); |
| | 1240 | |
| | 1241 | TIter next(trackSegments) ; |
| | 1242 | AliPHOSTrackSegment * ts ; |
| | 1243 | Int_t index = 0 ; |
| | 1244 | AliPHOSRecParticle * rp ; |
| | 1245 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
| | 1246 | // cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl; |
| | 1247 | new( (*recParticles)[index] ) AliPHOSRecParticle() ; |
| | 1248 | rp = (AliPHOSRecParticle *)recParticles->At(index) ; |
| | 1249 | rp->SetTrackSegment(index) ; |
| | 1250 | rp->SetIndexInList(index) ; |
| | 1251 | |
| | 1252 | AliPHOSEmcRecPoint * emc = 0 ; |
| | 1253 | if(ts->GetEmcIndex()>=0) |
| | 1254 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; |
| | 1255 | |
| | 1256 | AliPHOSCpvRecPoint * cpv = 0 ; |
| | 1257 | if(ts->GetCpvIndex()>=0) |
| | 1258 | cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
| | 1259 | |
| | 1260 | Int_t track = 0 ; |
| | 1261 | track = ts->GetTrackIndex() ; |
| | 1262 | |
| | 1263 | // Now set type (reconstructed) of the particle |
| | 1264 | |
| | 1265 | // Choose the cluster energy range |
| | 1266 | |
| | 1267 | if (!emc) { |
| | 1268 | AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; |
| | 1269 | } |
| | 1270 | |
| | 1271 | Float_t e = emc->GetEnergy() ; |
| | 1272 | |
| | 1273 | Float_t lambda[2] ; |
| | 1274 | emc->GetElipsAxis(lambda) ; |
| | 1275 | |
| | 1276 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ |
| | 1277 | // Looking PCA. Define and calculate the data (X), |
| | 1278 | // introduce in the function X2P that gives the components (P). |
| | 1279 | |
| | 1280 | Float_t spher = 0. ; |
| | 1281 | Float_t emaxdtotal = 0. ; |
| | 1282 | |
| | 1283 | if((lambda[0]+lambda[1])!=0) |
| | 1284 | spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); |
| | 1285 | |
| | 1286 | emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); |
| | 1287 | |
| | 1288 | fX[0] = lambda[0] ; |
| | 1289 | fX[1] = lambda[1] ; |
| | 1290 | fX[2] = emc->GetDispersion() ; |
| | 1291 | fX[3] = spher ; |
| | 1292 | fX[4] = emc->GetMultiplicity() ; |
| | 1293 | fX[5] = emaxdtotal ; |
| | 1294 | fX[6] = emc->GetCoreEnergy() ; |
| | 1295 | |
| | 1296 | fPrincipalPhoton->X2P(fX,fPPhoton); |
| | 1297 | fPrincipalPi0 ->X2P(fX,fPPi0); |
| | 1298 | |
| | 1299 | } |
| | 1300 | else{ |
| | 1301 | fPPhoton[0]=-100.0; //We do not accept clusters with |
| | 1302 | fPPhoton[1]=-100.0; //one cell as a photon-like |
| | 1303 | fPPi0[0] =-100.0; |
| | 1304 | fPPi0[1] =-100.0; |
| | 1305 | } |
| | 1306 | |
| | 1307 | Float_t time = emc->GetTime() ; |
| | 1308 | rp->SetTof(time) ; |
| | 1309 | |
| | 1310 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency |
| | 1311 | // are taken into account to set the particle identification) |
| | 1312 | for(Int_t effPur = 0; effPur < 3 ; effPur++){ |
| | 1313 | |
| | 1314 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, |
| | 1315 | // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) |
| | 1316 | // is set to 1 |
| | 1317 | if(GetCPVBit(emc, cpv, effPur,e) == 1 ){ |
| | 1318 | rp->SetPIDBit(effPur) ; |
| | 1319 | //cout<<"CPV bit "<<effPur<<endl; |
| | 1320 | } |
| | 1321 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th |
| | 1322 | // bit (depending on the efficiency-purity point )is set to 1 |
| | 1323 | if(time< (*fParameters)(3,effPur)) |
| | 1324 | rp->SetPIDBit(effPur+3) ; |
| | 1325 | |
| | 1326 | //Photon PCA |
| | 1327 | //If we are inside the ellipse, 7th, 8th or 9th |
| | 1328 | // bit (depending on the efficiency-purity point )is set to 1 |
| | 1329 | if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1) |
| | 1330 | rp->SetPIDBit(effPur+6) ; |
| | 1331 | |
| | 1332 | //Pi0 PCA |
| | 1333 | //If we are inside the ellipse, 10th, 11th or 12th |
| | 1334 | // bit (depending on the efficiency-purity point )is set to 1 |
| | 1335 | if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1) |
| | 1336 | rp->SetPIDBit(effPur+9) ; |
| | 1337 | } |
| | 1338 | if(GetHardPhotonBit(emc)) |
| | 1339 | rp->SetPIDBit(12) ; |
| | 1340 | if(GetHardPi0Bit (emc)) |
| | 1341 | rp->SetPIDBit(13) ; |
| | 1342 | |
| | 1343 | if(track >= 0) |
| | 1344 | rp->SetPIDBit(14) ; |
| | 1345 | |
| | 1346 | //Set momentum, energy and other parameters |
| | 1347 | Float_t encal = GetCalibratedEnergy(e); |
| | 1348 | TVector3 dir = GetMomentumDirection(emc,cpv) ; |
| | 1349 | dir.SetMag(encal) ; |
| | 1350 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; |
| | 1351 | rp->SetCalcMass(0); |
| | 1352 | rp->Name(); //If photon sets the particle pdg name to gamma |
| | 1353 | rp->SetProductionVertex(0,0,0,0); |
| | 1354 | rp->SetFirstMother(-1); |
| | 1355 | rp->SetLastMother(-1); |
| | 1356 | rp->SetFirstDaughter(-1); |
| | 1357 | rp->SetLastDaughter(-1); |
| | 1358 | rp->SetPolarisation(0,0,0); |
| | 1359 | //Set the position in global coordinate system from the RecPoint |
| | 1360 | AliPHOSGeometry * geom = gime->PHOSGeometry() ; |
| | 1361 | AliPHOSTrackSegment * ts = gime->TrackSegment(rp->GetPHOSTSIndex()) ; |
| | 1362 | AliPHOSEmcRecPoint * erp = gime->EmcRecPoint(ts->GetEmcIndex()) ; |
| | 1363 | TVector3 pos ; |
| | 1364 | geom->GetGlobal(erp, pos) ; |
| | 1365 | rp->SetPos(pos); |
| | 1366 | index++ ; |
| | 1367 | } |
| | 1368 | } |
| | 1369 | |
| | 1370 | //____________________________________________________________________________ |
| | 1371 | void AliPHOSPIDv1::Print(const Option_t *) const |
| | 1372 | { |
| | 1373 | // Print the parameters used for the particle type identification |
| | 1374 | |
| | 1375 | AliInfo("=============== AliPHOSPIDv1 ================") ; |
| | 1376 | printf("Making PID\n") ; |
| | 1377 | printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ; |
| | 1378 | printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ; |
| | 1379 | printf(" Matrix of Parameters: 14x4\n") ; |
| | 1380 | printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ; |
| | 1381 | printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ; |
| | 1382 | printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ; |
| | 1383 | printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ; |
| | 1384 | Printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; |
| | 1385 | fParameters->Print() ; |
| | 1386 | } |
| | 1387 | |
| | 1388 | |
| | 1389 | |
| | 1390 | //____________________________________________________________________________ |
| | 1391 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
| | 1392 | { |
| | 1393 | // Print table of reconstructed particles |
| | 1394 | |
| | 1395 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; |
| | 1396 | |
| | 1397 | TClonesArray * recParticles = gime->RecParticles() ; |
| | 1398 | |
| | 1399 | TString message ; |
| | 1400 | message = "\nevent " ; |
| | 1401 | message += gAlice->GetEvNumber() ; |
| | 1402 | message += " found " ; |
| | 1403 | message += recParticles->GetEntriesFast(); |
| | 1404 | message += " RecParticles\n" ; |
| | 1405 | |
| | 1406 | if(strstr(option,"all")) { // printing found TS |
| | 1407 | message += "\n PARTICLE Index \n" ; |
| | 1408 | |
| | 1409 | Int_t index ; |
| | 1410 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { |
| | 1411 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; |
| | 1412 | message += "\n" ; |
| | 1413 | message += rp->Name().Data() ; |
| | 1414 | message += " " ; |
| | 1415 | message += rp->GetIndexInList() ; |
| | 1416 | message += " " ; |
| | 1417 | message += rp->GetType() ; |
| | 1418 | } |
| | 1419 | } |
| | 1420 | AliInfo(message.Data() ) ; |
| | 1421 | } |
| | 1422 | |
| | 1423 | //____________________________________________________________________________ |
| | 1424 | void AliPHOSPIDv1::SetParameters() |
| | 1425 | { |
| | 1426 | // PCA : To do the Principal Components Analysis it is necessary |
| | 1427 | // the Principal file, which is opened here |
| | 1428 | fX = new double[7]; // Data for the PCA |
| | 1429 | fPPhoton = new double[7]; // Eigenvalues of the PCA |
| | 1430 | fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA |
| | 1431 | |
| | 1432 | // Read photon principals from the photon file |
| | 1433 | |
| | 1434 | fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; |
| | 1435 | TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ; |
| | 1436 | fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ; |
| | 1437 | f.Close() ; |
| | 1438 | |
| | 1439 | // Read pi0 principals from the pi0 file |
| | 1440 | |
| | 1441 | fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ; |
| | 1442 | TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ; |
| | 1443 | fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ; |
| | 1444 | fPi0.Close() ; |
| | 1445 | |
| | 1446 | // Open parameters file and initialization of the Parameters matrix. |
| | 1447 | // In the File Parameters.dat are all the parameters. These are introduced |
| | 1448 | // in a matrix of 16x4 |
| | 1449 | // |
| | 1450 | // All the parameters defined in this file are, in order of row: |
| | 1451 | // line 0 : calibration |
| | 1452 | // lines 1,2 : CPV rectangular cat for X and Z |
| | 1453 | // line 3 : TOF cut |
| | 1454 | // lines 4-8 : parameters to calculate photon PCA ellipse |
| | 1455 | // lines 9-13: parameters to calculate pi0 PCA ellipse |
| | 1456 | // lines 14-15: parameters to calculate border for high-pt photons and pi0 |
| | 1457 | |
| | 1458 | fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); |
| | 1459 | fParameters = new TMatrixF(16,4) ; |
| | 1460 | const Int_t kMaxLeng=255; |
| | 1461 | char string[kMaxLeng]; |
| | 1462 | |
| | 1463 | // Open a text file with PID parameters |
| | 1464 | FILE *fd = fopen(fFileNameParameters.Data(),"r"); |
| | 1465 | if (!fd) |
| | 1466 | AliFatal(Form("File %s with a PID parameters cannot be opened\n", |
| | 1467 | fFileNameParameters.Data())); |
| | 1468 | |
| | 1469 | Int_t i=0; |
| | 1470 | // Read parameter file line-by-line and skip empty line and comments |
| | 1471 | while (fgets(string,kMaxLeng,fd) != NULL) { |
| | 1472 | if (string[0] == '\n' ) continue; |
| | 1473 | if (string[0] == '!' ) continue; |
| | 1474 | sscanf(string, "%f %f %f %f", |
| | 1475 | &(*fParameters)(i,0), &(*fParameters)(i,1), |
| | 1476 | &(*fParameters)(i,2), &(*fParameters)(i,3)); |
| | 1477 | i++; |
| | 1478 | AliDebug(1, Form("SetParameters", "line %d: %s",i,string)); |
| | 1479 | } |
| | 1480 | fclose(fd); |
| | 1481 | } |
| | 1482 | |
| | 1483 | //____________________________________________________________________________ |
| | 1484 | void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) |
| | 1485 | { |
| | 1486 | // Set parameter "Calibration" i to a value param |
| | 1487 | if(i>2 || i<0) { |
| | 1488 | AliError(Form("Invalid parameter number: %d",i)); |
| | 1489 | } else |
| | 1490 | (*fParameters)(0,i) = param ; |
| | 1491 | } |
| | 1492 | |
| | 1493 | //____________________________________________________________________________ |
| | 1494 | void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) |
| | 1495 | { |
| | 1496 | // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on |
| | 1497 | // Purity-Efficiency point i |
| | 1498 | |
| | 1499 | if(i>2 || i<0) { |
| | 1500 | AliError(Form("Invalid parameter number: %d",i)); |
| | 1501 | } else { |
| | 1502 | axis.ToLower(); |
| | 1503 | if (axis == "x") (*fParameters)(1,i) = cut; |
| | 1504 | else if (axis == "z") (*fParameters)(2,i) = cut; |
| | 1505 | else { |
| | 1506 | AliError(Form("Invalid axis name: %s",axis.Data())); |
| | 1507 | } |
| | 1508 | } |
| | 1509 | } |
| | 1510 | |
| | 1511 | //____________________________________________________________________________ |
| | 1512 | void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) |
| | 1513 | { |
| | 1514 | // Set parameter "Hard photon boundary" i to a value param |
| | 1515 | if(i>4 || i<0) { |
| | 1516 | AliError(Form("Invalid parameter number: %d",i)); |
| | 1517 | } else |
| | 1518 | (*fParameters)(14,i) = param ; |
| | 1519 | } |
| | 1520 | |
| | 1521 | //____________________________________________________________________________ |
| | 1522 | void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) |
| | 1523 | { |
| | 1524 | // Set parameter "Hard pi0 boundary" i to a value param |
| | 1525 | if(i>1 || i<0) { |
| | 1526 | AliError(Form("Invalid parameter number: %d",i)); |
| | 1527 | } else |
| | 1528 | (*fParameters)(15,i) = param ; |
| | 1529 | } |
| | 1530 | |
| | 1531 | //_____________________________________________________________________________ |
| | 1532 | void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) |
| | 1533 | { |
| | 1534 | // Set the parameter TimeGate depending on Purity-Efficiency point i |
| | 1535 | if (i>2 || i<0) { |
| | 1536 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); |
| | 1537 | } else |
| | 1538 | (*fParameters)(3,i)= gate ; |
| | 1539 | } |
| | 1540 | |
| | 1541 | //_____________________________________________________________________________ |
| | 1542 | void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par) |
| | 1543 | { |
| | 1544 | // Set the parameter "i" that is needed to calculate the ellipse |
| | 1545 | // parameter "param" for a particle "particle" |
| | 1546 | |
| | 1547 | particle.ToLower(); |
| | 1548 | param. ToLower(); |
| | 1549 | Int_t p= -1; |
| | 1550 | Int_t offset=0; |
| | 1551 | |
| | 1552 | if (particle == "photon") offset=0; |
| | 1553 | else if (particle == "pi0") offset=5; |
| | 1554 | else |
| | 1555 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
| | 1556 | particle.Data())); |
| | 1557 | |
| | 1558 | if (param.Contains("a")) p=4+offset; |
| | 1559 | else if(param.Contains("b")) p=5+offset; |
| | 1560 | else if(param.Contains("c")) p=6+offset; |
| | 1561 | else if(param.Contains("x0"))p=7+offset; |
| | 1562 | else if(param.Contains("y0"))p=8+offset; |
| | 1563 | if((i>4)||(i<0)) { |
| | 1564 | AliError(Form("No parameter with index %d", i)) ; |
| | 1565 | } else if(p==-1) { |
| | 1566 | AliError(Form("No parameter with name %s", param.Data() )) ; |
| | 1567 | } else |
| | 1568 | (*fParameters)(p,i) = par ; |
| | 1569 | } |
| | 1570 | |
| | 1571 | //____________________________________________________________________________ |
| | 1572 | void AliPHOSPIDv1::Unload() |
| | 1573 | { |
| | 1574 | //Unloads RecPoints, Tracks and RecParticles |
| | 1575 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
| | 1576 | gime->PhosLoader()->UnloadRecPoints() ; |
| | 1577 | gime->PhosLoader()->UnloadTracks() ; |
| | 1578 | gime->PhosLoader()->UnloadRecParticles() ; |
| | 1579 | } |
| | 1580 | |
| | 1581 | //____________________________________________________________________________ |
| | 1582 | void AliPHOSPIDv1::WriteRecParticles() |
| | 1583 | { |
| | 1584 | //It writes reconstructed particles and pid to file |
| | 1585 | |
| | 1586 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; |
| | 1587 | |
| | 1588 | TClonesArray * recParticles = gime->RecParticles() ; |
| | 1589 | recParticles->Expand(recParticles->GetEntriesFast() ) ; |
| | 1590 | if(fWrite){ |
| | 1591 | TTree * treeP = gime->TreeP(); |
| | 1592 | |
| | 1593 | //First rp |
| | 1594 | Int_t bufferSize = 32000 ; |
| | 1595 | TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize); |
| | 1596 | rpBranch->SetTitle(BranchName()); |
| | 1597 | |
| | 1598 | rpBranch->Fill() ; |
| | 1599 | |
| | 1600 | gime->WriteRecParticles("OVERWRITE"); |
| | 1601 | gime->WritePID("OVERWRITE"); |
| | 1602 | } |
| | 1603 | } |
| | 1604 | |
| | 1605 | |
| | 1606 | //_______________________________________________________________________ |
| | 1607 | void AliPHOSPIDv1::SetInitPID(const Double_t *p) { |
| | 1608 | // Sets values for the initial population of each particle type |
| | 1609 | for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i]; |
| | 1610 | } |
| | 1611 | //_______________________________________________________________________ |
| | 1612 | void AliPHOSPIDv1::GetInitPID(Double_t *p) const { |
| | 1613 | // Gets values for the initial population of each particle type |
| | 1614 | for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i]; |
| | 1615 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental