From defc90400c9c69973b946d449d61f4a1918047d4 Mon Sep 17 00:00:00 2001 From: cblume Date: Mon, 2 Jan 2006 07:51:48 +0000 Subject: [PATCH] New offline PID code by Prashant --- TRD/AliTRDpidESD.cxx | 10 +- TRD/AliTRDpidESD.h | 5 +- TRD/AliTRDprobdist.cxx | 3427 +++------------------------------- TRD/AliTRDprobdist.h | 98 +- TRD/ConfigPID.C | 797 ++++++++ TRD/TRDdEdxHistogramsV1.root | Bin 0 -> 53434 bytes TRD/histTRDpid.C | 315 ++++ 7 files changed, 1386 insertions(+), 3266 deletions(-) create mode 100644 TRD/ConfigPID.C create mode 100644 TRD/TRDdEdxHistogramsV1.root create mode 100644 TRD/histTRDpid.C diff --git a/TRD/AliTRDpidESD.cxx b/TRD/AliTRDpidESD.cxx index 8f9e4e16ec6..ac899f27f0c 100644 --- a/TRD/AliTRDpidESD.cxx +++ b/TRD/AliTRDpidESD.cxx @@ -93,8 +93,14 @@ Int_t AliTRDpidESD::MakePID(AliESD *event) // // The class AliTRDprobdist contains precalculated prob dis. AliTRDprobdist *pd = new AliTRDprobdist(); - pd->SetADCNorm(1.0); - // pd->SetADCNorm(.72); + pd->SetADCNorm(1.0); // The factor is the ratio of Mean of pi charge dist. + // for the New TRD code divided by the Mean of pi charge + // dist. given in AliTRDprobdist object + + // Example to get mean for particle 2 (pi) and momentum number 4 (2 GeV) + // printf("%.2f \n", pd->GetMean(2, 4)); + // Example of use of Copy Constructor + // AliTRDprobdist *pd1 = new AliTRDprobdist(*pd); Int_t ntrk=event->GetNumberOfTracks(); for (Int_t i=0; i @@ -16,7 +15,7 @@ class AliTRDpidESD { public: AliTRDpidESD(Double_t *param); virtual ~AliTRDpidESD() {} - Int_t MakePID(AliESD *event); + static Int_t MakePID(AliESD *event); static Double_t Bethe(Double_t bg); private: Double_t fMIP; // dEdx for MIP diff --git a/TRD/AliTRDprobdist.cxx b/TRD/AliTRDprobdist.cxx index 58c6d21ebcb..7de83cf7a85 100644 --- a/TRD/AliTRDprobdist.cxx +++ b/TRD/AliTRDprobdist.cxx @@ -26,194 +26,294 @@ ClassImp(AliTRDprobdist) //_________________________________________________________________________ -AliTRDprobdist::AliTRDprobdist(Int_t multiplicity) +AliTRDprobdist::AliTRDprobdist(Char_t *responseFile) { // // The main constructor // - if (multiplicity == 1) FillData(); - // if (multiplicity == 2000) FillData2000(); - // if (multiplicity == 4000) FillData4000(); - // if (multiplicity == 6000) FillData6000(); - // if (multiplicity == 8000) FillData8000(); + Char_t *partName[5] = {"electron", "muon", "pion", "kaon", "proton"}; + for(Int_t i=0; i<5; i++) { + fpartName[i] = partName[i]; + } + // ADC Gain normalization + fADCNorm=1.0; + fNMom=kNMom; + // Set track momenta for which response functions are available + Double_t trackMomentum[kNMom]= {0.6, 0.8, 1, 1.5, 2, 3, 4, 5, 6, 8, 10}; + for(Int_t imom=0; imomIsOpen()) { + Error("AliTRDprobdist", "opening TRD histgram file %s failed", responseFile); + return kFALSE; } - return integral/norm; + + // Read histograms + Char_t text[200]; + for (Int_t imom = 0; imom < kNMom; imom++) { + sprintf(text,"h1dEdxEL%01d",imom+1); + fh1dEdxEL[imom] = (TH1F*)histFile->Get(text); + fh1dEdxEL[imom]->Scale(1.0/fh1dEdxEL[imom]->Integral()); + + sprintf(text,"h1dEdxPI%01d",imom+1); + fh1dEdxPI[imom] = (TH1F*)histFile->Get(text); + fh1dEdxPI[imom]->Scale(1.0/fh1dEdxPI[imom]->Integral()); + + sprintf(text,"h1dEdxMU%01d",imom+1); + fh1dEdxMU[imom] = (TH1F*)histFile->Get(text); + fh1dEdxMU[imom]->Scale(1.0/fh1dEdxMU[imom]->Integral()); + + sprintf(text,"h1dEdxKA%01d",imom+1); + fh1dEdxKA[imom] = (TH1F*)histFile->Get(text); + fh1dEdxKA[imom]->Scale(1.0/fh1dEdxKA[imom]->Integral()); + + sprintf(text,"h1dEdxPR%01d",imom+1); + fh1dEdxPR[imom] = (TH1F*)histFile->Get(text); + fh1dEdxPR[imom]->Scale(1.0/fh1dEdxPR[imom]->Integral()); + + sprintf(text,"h1MaxTimBinEL%01d",imom+1); + fh1MaxTimBinEL[imom] = (TH1F*)histFile->Get(text); + fh1MaxTimBinEL[imom]->Scale(1.0/fh1MaxTimBinEL[imom]->Integral()); + + sprintf(text,"h1MaxTimBinPI%01d",imom+1); + fh1MaxTimBinPI[imom] = (TH1F*)histFile->Get(text); + fh1MaxTimBinPI[imom]->Scale(1.0/fh1MaxTimBinPI[imom]->Integral()); + } + // Number of bins and bin size + fNbins = fh1dEdxPI[1]->GetNbinsX(); + fBinSize = fh1dEdxPI[1]->GetBinWidth(1); + return kTRUE; } -//_________________________________________________________________________ -Double_t AliTRDprobdist::GetMeanEL(Int_t ip) const +AliTRDprobdist::AliTRDprobdist(const AliTRDprobdist& pd):TNamed() { // - // Gets mean of de/dx dist. of e - Double_t integral=0.; - Double_t norm=0.; - for(Int_t ie=0; ieGetMean(); + if(k==1) return fh1dEdxMU[ip]->GetMean(); + if(k==2) return fh1dEdxPI[ip]->GetMean(); + if(k==3) return fh1dEdxKA[ip]->GetMean(); + if(k==4) return fh1dEdxPR[ip]->GetMean(); + return fh1dEdxPR[ip]->GetMean(); } //_________________________________________________________________________ -Double_t AliTRDprobdist::GetNormalizationEL(Int_t ip) const +Double_t AliTRDprobdist::GetNormalization(Int_t k, Int_t ip) const { // // Gets Normalization of de/dx dist. of e - Double_t integral=0.; - for(Int_t ie=0; ieIntegral(); + if(k==1) return fh1dEdxMU[ip]->Integral(); + if(k==2) return fh1dEdxPI[ip]->Integral(); + if(k==3) return fh1dEdxKA[ip]->Integral(); + if(k==4) return fh1dEdxPR[ip]->Integral(); + return fh1dEdxPR[ip]->Integral(); +} + +TH1F* AliTRDprobdist::GetHistogram(Int_t k, Int_t ip) const +{ + // + // + printf("Histogram for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]); + if(k==0) return fh1dEdxEL[ip]; + if(k==1) return fh1dEdxMU[ip]; + if(k==2) return fh1dEdxPI[ip]; + if(k==3) return fh1dEdxKA[ip]; + if(k==4) return fh1dEdxPR[ip]; + return fh1dEdxPR[ip]; } //_________________________________________________________________________ -Double_t AliTRDprobdist::GetProbability(Int_t k, Double_t mom, Double_t dedx) const +Double_t AliTRDprobdist::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const { // // Gets the Probability of having dedx at a given momentum (mom) // and particle type k (0 for e) and (2 for pi) // from the precalculated de/dx distributions Double_t probability = 1.0; - Int_t iEnBin= ((Int_t) (dedx/fEnBinSize)); - if(iEnBin > fNEbins-1) iEnBin = fNEbins-1; + Double_t dedx = dedx1/fADCNorm; + Int_t iEnBin= ((Int_t) (dedx/fBinSize+1)); + if(iEnBin > fNbins) iEnBin = fNbins; - if(k==0){ // electron + ////Electron////////////////////////// + if(k==0){ // electron Double_t slop; + // Lower limit if(mom<=fTrackMomentum[0]) { - slop=(fProbEL[1][iEnBin]-fProbEL[0][iEnBin])/(fTrackMomentum[1] - fTrackMomentum[0]); - probability= fProbEL[0][iEnBin] + slop*(mom-fTrackMomentum[0]); + slop=(fh1dEdxEL[1]->GetBinContent(iEnBin)-fh1dEdxEL[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]); + probability= fh1dEdxEL[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]); return probability; } + // Upper Limit if(mom>=fTrackMomentum[fNMom-1]) { - slop=(fProbEL[fNMom-1][iEnBin]-fProbEL[fNMom-2][iEnBin])/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); - probability= fProbEL[fNMom-1][iEnBin] + slop*(mom-fTrackMomentum[fNMom-1]); + slop=(fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxEL[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); + probability= fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]); return probability; } + // In the range + for(Int_t ip=1; ipGetBinContent(iEnBin)-fh1dEdxEL[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]); + // Linear Interpolation + probability= fh1dEdxEL[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]); + return probability; + } } - if(k==2){ // pion + ////Pion////////////////////////// + if(k==2){ // Pion Double_t slop; + // Lower limit if(mom<=fTrackMomentum[0]) { - slop=(fProbPI[1][iEnBin]-fProbPI[0][iEnBin])/(fTrackMomentum[1] - fTrackMomentum[0]); - probability= fProbPI[0][iEnBin] + slop*(mom-fTrackMomentum[0]); + slop=(fh1dEdxPI[1]->GetBinContent(iEnBin)-fh1dEdxPI[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]); + probability= fh1dEdxPI[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]); return probability; } + // Upper Limit if(mom>=fTrackMomentum[fNMom-1]) { - slop=(fProbPI[fNMom-1][iEnBin]-fProbPI[fNMom-2][iEnBin])/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); - probability= fProbPI[fNMom-1][iEnBin] + slop*(mom-fTrackMomentum[fNMom-1]); + slop=(fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPI[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); + probability= fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]); return probability; } + // In the range + for(Int_t ip=1; ipGetBinContent(iEnBin)-fh1dEdxPI[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]); + // Linear Interpolation + probability= fh1dEdxPI[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]); + return probability; + } } - if(k==1){ // muon + ////Muon////////////////////////// + if(k==1){ // Muon Double_t slop; + // Lower limit if(mom<=fTrackMomentum[0]) { - slop=(fProbMU[1][iEnBin]-fProbMU[0][iEnBin])/(fTrackMomentum[1] - fTrackMomentum[0]); - probability= fProbMU[0][iEnBin] + slop*(mom-fTrackMomentum[0]); + slop=(fh1dEdxMU[1]->GetBinContent(iEnBin)-fh1dEdxMU[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]); + probability= fh1dEdxMU[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]); return probability; } + // Upper Limit if(mom>=fTrackMomentum[fNMom-1]) { - slop=(fProbMU[fNMom-1][iEnBin]-fProbMU[fNMom-2][iEnBin])/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); - probability= fProbMU[fNMom-1][iEnBin] + slop*(mom-fTrackMomentum[fNMom-1]); + slop=(fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxMU[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); + probability= fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]); return probability; } + // In the range + for(Int_t ip=1; ipGetBinContent(iEnBin)-fh1dEdxMU[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]); + // Linear Interpolation + probability= fh1dEdxMU[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]); + return probability; + } } - if(k==3){ // kaon + ////Kaon////////////////////////// + if(k==3){ // Kaon Double_t slop; + // Lower limit if(mom<=fTrackMomentum[0]) { - slop=(fProbKA[1][iEnBin]-fProbKA[0][iEnBin])/(fTrackMomentum[1] - fTrackMomentum[0]); - probability= fProbKA[0][iEnBin] + slop*(mom-fTrackMomentum[0]); + slop=(fh1dEdxKA[1]->GetBinContent(iEnBin)-fh1dEdxKA[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]); + probability= fh1dEdxKA[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]); return probability; } + // Upper Limit if(mom>=fTrackMomentum[fNMom-1]) { - slop=(fProbKA[fNMom-1][iEnBin]-fProbKA[fNMom-2][iEnBin])/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); - probability= fProbKA[fNMom-1][iEnBin] + slop*(mom-fTrackMomentum[fNMom-1]); + slop=(fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxKA[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); + probability= fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]); return probability; } + // In the range + for(Int_t ip=1; ipGetBinContent(iEnBin)-fh1dEdxKA[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]); + // Linear Interpolation + probability= fh1dEdxKA[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]); + return probability; + } } - if(k==4){ // proton + ////Proton////////////////////////// + if(k==4){ // Proton Double_t slop; + // Lower limit if(mom<=fTrackMomentum[0]) { - slop=(fProbPR[1][iEnBin]-fProbPR[0][iEnBin])/(fTrackMomentum[1] - fTrackMomentum[0]); - probability= fProbPR[0][iEnBin] + slop*(mom-fTrackMomentum[0]); + slop=(fh1dEdxPR[1]->GetBinContent(iEnBin)-fh1dEdxPR[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]); + probability= fh1dEdxPR[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]); return probability; } + // Upper Limit if(mom>=fTrackMomentum[fNMom-1]) { - slop=(fProbPR[fNMom-1][iEnBin]-fProbPR[fNMom-2][iEnBin])/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); - probability= fProbPR[fNMom-1][iEnBin] + slop*(mom-fTrackMomentum[fNMom-1]); + slop=(fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPR[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]); + probability= fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]); return probability; } + // In the range + for(Int_t ip=1; ipGetBinContent(iEnBin)-fh1dEdxPR[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]); + // Linear Interpolation + probability= fh1dEdxPR[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]); + return probability; + } } - if(k==0) // electron - for(Int_t ip=1; ip=fTrackMomentum[fNMom-1]) probabilityT = fProbELT[fNMom-1][iTBin]; + if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinEL[0]->GetBinContent(iTBin); + if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinEL[fNMom-1]->GetBinContent(iTBin); } if(k==1||k==2||k==3||k==4){ // pion - if(mom<=fTrackMomentum[0]) probabilityT = fProbPIT[0][iTBin]; - if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fProbPIT[fNMom-1][iTBin]; + if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinPI[0]->GetBinContent(iTBin); + if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinPI[fNMom-1]->GetBinContent(iTBin); } if(k==0) // electron for(Int_t ip=1; ipGetBinContent(iTBin)-fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]); // Linear Interpolation - probabilityT= fProbELT[ip-1][iTBin] + slop*(mom-fTrackMomentum[ip-1]); + probabilityT= fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]); return probabilityT; } if(k==1||k==2||k==3||k==4) // pion and other particles for(Int_t ip=1; ipGetBinContent(iTBin)-fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]); // Linear Interpolation - probabilityT= fProbPIT[ip-1][iTBin] + slop*(mom-fTrackMomentum[ip-1]); + probabilityT= fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]); return probabilityT; } return probabilityT; } -void AliTRDprobdist::SetADCNorm(Double_t norm) -{ - // - // Setting the multiplicative constant to Energy loss scale of ADC - fADCNorm=norm; - fEnBinSize=norm*fEnBinSize; - for(Int_t ie=0; ie -----------------------------------------------------------------*/ -#include "TObject.h" +#include +#include +#include -const Int_t kNo_Mom=11; -const Int_t kNo_EnBins=250; -const Int_t kNo_TBins=20; +const Int_t kNMom=11; -class AliTRDprobdist : public TObject { +class AliTRDprobdist : public TNamed { public: - AliTRDprobdist(Int_t multiplicity=1); // multiplicity can take - // values 1, 2000, 4000, 6000, 8000 - // ~AliTRDprobdist(); - void FillData(); // Fills Data - // void FillData2000(); // Fills Data of multiplicity 2000 - // void FillData4000(); // Fills Data of multiplicity 4000 - // void FillData6000(); // Fills Data of multiplicity 6000 - // void FillData8000(); // Fills Data of multiplicity 8000 - Double_t GetBM(Int_t ip) const { return fTrackMomentum[ip];} - Double_t GetMeanPI(Int_t ip) const; // Gets mean of de/dx dist. of pi - Double_t GetMeanEL(Int_t ip) const; // Gets mean of de/dx dist. of e - Double_t GetNormalizationPI(Int_t ip) const; // Gets Norm. of de/dx dist. of pi - Double_t GetNormalizationEL(Int_t ip) const; // Gets Norm. of de/dx dist. of e - Double_t GetProbability(Int_t k, Double_t mom, Double_t dedx) const; - // Gets the Probability of having dedx - Double_t GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const; - // Gets the Probability of having timbin - void SetADCNorm(Double_t norm); - void GetData(Int_t ip, Double_t *ebin, Double_t *ppi, Double_t *pel) const { - for(Int_t ie=0; ieLoad("libgeant321"); +// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\ +// -I$ALICE_ROOT -I$ALICE/geant3/TGeant3"); +// root [0] .x grun.C(1,"ConfigPPR.C++") +// Prashant Shukla (shukla@physi.uni-heidelberg.de) + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include "STEER/AliRunLoader.h" +#include "STEER/AliRun.h" +#include "STEER/AliConfig.h" +#include "STEER/AliGenerator.h" +#include "STEER/AliLog.h" +#include "PYTHIA6/AliDecayerPythia.h" +#include "EVGEN/AliGenHIJINGpara.h" +#include "THijing/AliGenHijing.h" +#include "EVGEN/AliGenCocktail.h" +#include "EVGEN/AliGenSlowNucleons.h" +#include "EVGEN/AliSlowNucleonModelExp.h" +#include "EVGEN/AliGenParam.h" +#include "EVGEN/AliGenMUONlib.h" +#include "EVGEN/AliGenSTRANGElib.h" +#include "EVGEN/AliGenMUONCocktail.h" +#include "EVGEN/AliGenCocktail.h" +#include "EVGEN/AliGenGeVSim.h" +#include "EVGEN/AliGeVSimParticle.h" +#include "PYTHIA6/AliGenPythia.h" +#include "STEER/AliMagFMaps.h" +#include "STRUCT/AliBODY.h" +#include "STRUCT/AliMAG.h" +#include "STRUCT/AliABSOv0.h" +#include "STRUCT/AliDIPOv2.h" +#include "STRUCT/AliHALL.h" +#include "STRUCT/AliFRAMEv2.h" +#include "STRUCT/AliSHILv2.h" +#include "STRUCT/AliPIPEv0.h" +#include "ITS/AliITSvPPRasymmFMD.h" +#include "TPC/AliTPCv2.h" +#include "TOF/AliTOFv4T0.h" +#include "RICH/AliRICHv1.h" +#include "ZDC/AliZDCv2.h" +#include "TRD/AliTRDv1.h" +#include "FMD/AliFMDv1.h" +#include "MUON/AliMUONv1.h" +#include "MUON/AliMUONSt1GeometryBuilderV2.h" +#include "MUON/AliMUONSt2GeometryBuilderV2.h" +#include "MUON/AliMUONSlatGeometryBuilder.h" +#include "MUON/AliMUONTriggerGeometryBuilder.h" +#include "PHOS/AliPHOSv1.h" +#include "PMD/AliPMDv1.h" +#include "START/AliSTARTv1.h" +#include "EMCAL/AliEMCALv1.h" +#include "CRT/AliCRTv0.h" +#include "VZERO/AliVZEROv5.h" +#endif + +enum PprRun_t +{ + test50, + kParam_8000, kParam_4000, kParam_2000, + kHijing_cent1, kHijing_cent2, + kHijing_per1, kHijing_per2, kHijing_per3, kHijing_per4, kHijing_per5, + kCocktailTRD, kpieTRD +}; + +const char* pprRunName[] = { + "test50", + "kParam_8000", "kParam_4000", "kParam_2000", + "kHijing_cent1", "kHijing_cent2", + "kHijing_per1", "kHijing_per2", "kHijing_per3", "kHijing_per4", + "kHijing_per5", + "kCocktailTRD", "kpieTRD" +}; + +enum PprGeo_t +{ + kHoles, kNoHoles +}; + +enum PprRad_t +{ + kGluonRadiation, kNoGluonRadiation +}; + +enum PprMag_t +{ + k2kG, k4kG, k5kG +}; + + +// This part for configuration +//static PprRun_t srun = test50; +static PprRun_t srun = kpieTRD; +static PprGeo_t sgeo = kNoHoles; +static PprRad_t srad = kGluonRadiation; +static PprMag_t smag = k5kG; +static Int_t sseed = 0; //Set 0 to use the current time + +// Comment line +static TString comment; + +// Functions +Float_t EtaToTheta(Float_t arg); +AliGenerator* GeneratorFactory(PprRun_t srun); +AliGenHijing* HijingStandard(); +void ProcessEnvironmentVars(); + +void Config() +{ + // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00 + // Theta range given through pseudorapidity limits 22/6/2001 + + // Get settings from environment variables + ProcessEnvironmentVars(); + + // Set Random Number seed + gRandom->SetSeed(sseed); + cout<<"Seed for random number generation= "<GetSeed()<Load("libgeant321"); +#endif + + new TGeant3TGeo("C++ Interface to Geant3"); + + AliRunLoader* rl=0x0; + + cout<<"Config.C: Creating Run Loader ..."<Fatal("Config.C","Can not instatiate the Run Loader"); + return; + } + rl->SetCompressionLevel(2); + rl->SetNumberOfEventsPerFile(3); + gAlice->SetRunLoader(rl); + + // + // Set External decayer + AliDecayer *decayer = new AliDecayerPythia(); + + decayer->SetForceDecay(kAll); + + decayer->Init(); + gMC->SetExternalDecayer(decayer); + // + // + //======================================================================= + // + //======================================================================= + // ************* STEERING parameters FOR ALICE SIMULATION ************** + // --- Specify event type to be tracked through the ALICE setup + // --- All positions are in cm, angles in degrees, and P and E in GeV + + gMC->SetProcess("DCAY",1); + gMC->SetProcess("PAIR",1); + gMC->SetProcess("COMP",1); + gMC->SetProcess("PHOT",1); + gMC->SetProcess("PFIS",0); + gMC->SetProcess("DRAY",0); + gMC->SetProcess("ANNI",1); + gMC->SetProcess("BREM",1); + gMC->SetProcess("MUNU",1); + gMC->SetProcess("CKOV",1); + gMC->SetProcess("HADR",1); + gMC->SetProcess("LOSS",2); + gMC->SetProcess("MULS",1); + gMC->SetProcess("RAYL",1); + + Float_t cut = 1.e-3; // 1MeV cut by default + Float_t tofmax = 1.e10; + + gMC->SetCut("CUTGAM", cut); + gMC->SetCut("CUTELE", cut); + gMC->SetCut("CUTNEU", cut); + gMC->SetCut("CUTHAD", cut); + gMC->SetCut("CUTMUO", cut); + gMC->SetCut("BCUTE", cut); + gMC->SetCut("BCUTM", cut); + gMC->SetCut("DCUTE", cut); + gMC->SetCut("DCUTM", cut); + gMC->SetCut("PPCUTM", cut); + gMC->SetCut("TOFMAX", tofmax); + + // Debug and log level + // AliLog::SetGlobalDebugLevel(0); + // AliLog::SetGlobalLogLevel(AliLog::kError); + + // Generator Configuration + AliGenerator* gener = GeneratorFactory(srun); + gener->SetOrigin(0, 0, 0); // vertex position + gener->SetSigma(0, 0, 5.3); // Sigma in (X,Y,Z) (cm) on IP position + gener->SetCutVertexZ(1.); // Truncate at 1 sigma + gener->SetVertexSmear(kPerEvent); + gener->SetTrackingFlag(1); + gener->Init(); + + if (smag == k2kG) { + comment = comment.Append(" | L3 field 0.2 T"); + } else if (smag == k4kG) { + comment = comment.Append(" | L3 field 0.4 T"); + } else if (smag == k5kG) { + comment = comment.Append(" | L3 field 0.5 T"); + } + + + if (srad == kGluonRadiation) + { + comment = comment.Append(" | Gluon Radiation On"); + + } else { + comment = comment.Append(" | Gluon Radiation Off"); + } + + if (sgeo == kHoles) + { + comment = comment.Append(" | Holes for PHOS/RICH"); + + } else { + comment = comment.Append(" | No holes for PHOS/RICH"); + } + + printf("\n \n Comment: %s \n \n", comment.Data()); + + +// Field (L3 0.4 T) + AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., smag); + // AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 0, 10., smag); // B = 0 + field->SetL3ConstField(0); //Using const. field in the barrel + rl->CdGAFile(); + gAlice->SetField(field); +// + Int_t iABSO = 1; + Int_t iDIPO = 1; + Int_t iFMD = 1; + Int_t iFRAME = 1; + Int_t iHALL = 1; + Int_t iITS = 1; + Int_t iMAG = 1; + Int_t iMUON = 0; + Int_t iPHOS = 1; + Int_t iPIPE = 1; + Int_t iPMD = 0; + Int_t iRICH = 0; + Int_t iSHIL = 1; + Int_t iSTART = 1; + Int_t iTOF = 0; + Int_t iTPC = 1; + Int_t iTRD = 1; + Int_t iZDC = 1; + Int_t iEMCAL = 0; + Int_t iVZERO = 1; + Int_t iCRT = 0; + + //=================== Alice BODY parameters ============================= + AliBODY *BODY = new AliBODY("BODY", "Alice envelop"); + + + if (iMAG) + { + //=================== MAG parameters ============================ + // --- Start with Magnet since detector layouts may be depending --- + // --- on the selected Magnet dimensions --- + AliMAG *MAG = new AliMAG("MAG", "Magnet"); + } + + + if (iABSO) + { + //=================== ABSO parameters ============================ + AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber"); + } + + if (iDIPO) + { + //=================== DIPO parameters ============================ + + AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2"); + } + + if (iHALL) + { + //=================== HALL parameters ============================ + + AliHALL *HALL = new AliHALL("HALL", "Alice Hall"); + } + + + if (iFRAME) + { + //=================== FRAME parameters ============================ + + AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame"); + if (sgeo == kHoles) { + FRAME->SetHoles(1); + } else { + FRAME->SetHoles(0); + } + } + + if (iSHIL) + { + //=================== SHIL parameters ============================ + + AliSHIL *SHIL = new AliSHILv2("SHIL", "Shielding Version 2"); + } + + + if (iPIPE) + { + //=================== PIPE parameters ============================ + + AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe"); + } + + if(iITS) { + + //=================== ITS parameters ============================ + // + // As the innermost detector in ALICE, the Inner Tracking System "impacts" on + // almost all other detectors. This involves the fact that the ITS geometry + // still has several options to be followed in parallel in order to determine + // the best set-up which minimizes the induced background. All the geometries + // available to date are described in the following. Read carefully the comments + // and use the default version (the only one uncommented) unless you are making + // comparisons and you know what you are doing. In this case just uncomment the + // ITS geometry you want to use and run Aliroot. + // + // Detailed geometries: + // + // + //AliITS *ITS = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services"); + // + //AliITS *ITS = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services"); + // + AliITSvPPRasymmFMD *ITS = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services"); + ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer + ITS->SetReadDet(kTRUE); // don't touch this parameter if you're not an ITS developer + // ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det"); // don't touch this parameter if you're not an ITS developer + ITS->SetThicknessDet1(200.); // detector thickness on layer 1 must be in the range [100,300] + ITS->SetThicknessDet2(200.); // detector thickness on layer 2 must be in the range [100,300] + ITS->SetThicknessChip1(200.); // chip thickness on layer 1 must be in the range [150,300] + ITS->SetThicknessChip2(200.); // chip thickness on layer 2 must be in the range [150,300] + ITS->SetRails(0); // 1 --> rails in ; 0 --> rails out + ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon + + // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful + // for reconstruction !): + // + // + //AliITSvPPRcoarseasymm *ITS = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services"); + //ITS->SetRails(0); // 1 --> rails in ; 0 --> rails out + //ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon + // + //AliITS *ITS = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services"); + //ITS->SetRails(0); // 1 --> rails in ; 0 --> rails out + //ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon + // + // + // + // Geant3 <-> EUCLID conversion + // ============================ + // + // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and + // media to two ASCII files (called by default ITSgeometry.euc and + // ITSgeometry.tme) in a format understandable to the CAD system EUCLID. + // The default (=0) means that you dont want to use this facility. + // + ITS->SetEUCLID(0); + } + + if (iTPC) + { + //============================ TPC parameters ================================ + // --- This allows the user to specify sectors for the SLOW (TPC geometry 2) + // --- Simulator. SecAL (SecAU) <0 means that ALL lower (upper) + // --- sectors are specified, any value other than that requires at least one + // --- sector (lower or upper)to be specified! + // --- Reminder: sectors 1-24 are lower sectors (1-12 -> z>0, 13-24 -> z<0) + // --- sectors 25-72 are the upper ones (25-48 -> z>0, 49-72 -> z<0) + // --- SecLows - number of lower sectors specified (up to 6) + // --- SecUps - number of upper sectors specified (up to 12) + // --- Sens - sensitive strips for the Slow Simulator !!! + // --- This does NOT work if all S or L-sectors are specified, i.e. + // --- if SecAL or SecAU < 0 + // + // + //----------------------------------------------------------------------------- + + // gROOT->LoadMacro("SetTPCParam.C"); + // AliTPCParam *param = SetTPCParam(); + AliTPC *TPC = new AliTPCv2("TPC", "Default"); + + // All sectors included + TPC->SetSecAL(-1); + TPC->SetSecAU(-1); + + } + + + if (iTOF) { + //=================== TOF parameters ============================ + AliTOF *TOF = new AliTOFv4T0("TOF", "normal TOF"); + } + + + if (iRICH) + { + //=================== RICH parameters =========================== + AliRICH *RICH = new AliRICHv1("RICH", "normal RICH"); + + } + + + if (iZDC) + { + //=================== ZDC parameters ============================ + + AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC"); + } + + if (iTRD) + { + //=================== TRD parameters ============================ + + AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator"); + + // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 90% Xe + 10% CO2) + TRD->SetGasMix(1); + if (sgeo == kHoles) { + // With hole in front of PHOS + TRD->SetPHOShole(); + // With hole in front of RICH + TRD->SetRICHhole(); + } + // Switch on TR + AliTRDsim *TRDsim = TRD->CreateTR(); + } + + if (iFMD) + { + //=================== FMD parameters ============================ + AliFMD *FMD = new AliFMDv1("FMD", "normal FMD"); + } + + if (iMUON) + { + //=================== MUON parameters =========================== + // New MUONv1 version (geometry defined via builders) + AliMUON *MUON = new AliMUONv1("MUON", "default"); + ((AliMUONv1*)MUON)->SetStepManagerVersionDE(true); + MUON->AddGeometryBuilder(new AliMUONSt1GeometryBuilderV2(MUON)); + MUON->AddGeometryBuilder(new AliMUONSt2GeometryBuilderV2(MUON)); + MUON->AddGeometryBuilder(new AliMUONSlatGeometryBuilder(MUON)); + MUON->AddGeometryBuilder(new AliMUONTriggerGeometryBuilder(MUON)); + } + //=================== PHOS parameters =========================== + + if (iPHOS) + { + AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP"); + } + + + if (iPMD) + { + //=================== PMD parameters ============================ + AliPMD *PMD = new AliPMDv1("PMD", "normal PMD"); + } + + if (iSTART) + { + //=================== START parameters ============================ + AliSTART *START = new AliSTARTv1("START", "START Detector"); + } + + if (iEMCAL) + { + //=================== EMCAL parameters ============================ + AliEMCAL *EMCAL = new AliEMCALv1("EMCAL", "EMCAL_55_25"); + } + + if (iCRT) + { + //=================== CRT parameters ============================ + AliCRT *CRT = new AliCRTv0("CRT", "normal ACORDE"); + } + + if (iVZERO) + { + //=================== CRT parameters ============================ + AliVZERO *VZERO = new AliVZEROv5("VZERO", "normal VZERO"); + } + + +} + +Float_t EtaToTheta(Float_t arg){ + return (180./TMath::Pi())*2.*atan(exp(-arg)); +} + + + +AliGenerator* GeneratorFactory(PprRun_t srun) { + Int_t isw = 3; + if (srad == kNoGluonRadiation) isw = 0; + + + AliGenerator * gGener = 0x0; + switch (srun) { + case test50: + { + comment = comment.Append(":HIJINGparam test 50 particles"); + AliGenHIJINGpara *gener = new AliGenHIJINGpara(50); + gener->SetMomentumRange(0, 999999.); + gener->SetPhiRange(0., 360.); + // Set pseudorapidity range from -8 to 8. + Float_t thmin = EtaToTheta(8); // theta min. <---> eta max + Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min + gener->SetThetaRange(thmin,thmax); + gGener=gener; + } + break; + case kParam_8000: + { + comment = comment.Append(":HIJINGparam N=8000"); + AliGenHIJINGpara *gener = new AliGenHIJINGpara(86030); + gener->SetMomentumRange(0, 999999.); + gener->SetPhiRange(0., 360.); + // Set pseudorapidity range from -8 to 8. + Float_t thmin = EtaToTheta(8); // theta min. <---> eta max + Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min + gener->SetThetaRange(thmin,thmax); + gGener=gener; + } + break; + case kParam_4000: + { + comment = comment.Append("HIJINGparam N=4000"); + AliGenHIJINGpara *gener = new AliGenHIJINGpara(43015); + gener->SetMomentumRange(0, 999999.); + gener->SetPhiRange(0., 360.); + // Set pseudorapidity range from -8 to 8. + Float_t thmin = EtaToTheta(8); // theta min. <---> eta max + Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min + gener->SetThetaRange(thmin,thmax); + gGener=gener; + } + break; + case kParam_2000: + { + comment = comment.Append("HIJINGparam N=2000"); + AliGenHIJINGpara *gener = new AliGenHIJINGpara(21507); + gener->SetMomentumRange(0, 999999.); + gener->SetPhiRange(0., 360.); + // Set pseudorapidity range from -8 to 8. + Float_t thmin = EtaToTheta(8); // theta min. <---> eta max + Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min + gener->SetThetaRange(thmin,thmax); + gGener=gener; + } + break; +// +// Hijing Central +// + case kHijing_cent1: + { + comment = comment.Append("HIJING cent1"); + AliGenHijing *gener = HijingStandard(); +// impact parameter range + gener->SetImpactParameterRange(0., 5.); + gGener=gener; + } + break; + case kHijing_cent2: + { + comment = comment.Append("HIJING cent2"); + AliGenHijing *gener = HijingStandard(); +// impact parameter range + gener->SetImpactParameterRange(0., 2.); + gGener=gener; + } + break; +// +// Hijing Peripheral +// + case kHijing_per1: + { + comment = comment.Append("HIJING per1"); + AliGenHijing *gener = HijingStandard(); +// impact parameter range + gener->SetImpactParameterRange(5., 8.6); + gGener=gener; + } + break; + case kHijing_per2: + { + comment = comment.Append("HIJING per2"); + AliGenHijing *gener = HijingStandard(); +// impact parameter range + gener->SetImpactParameterRange(8.6, 11.2); + gGener=gener; + } + break; + case kHijing_per3: + { + comment = comment.Append("HIJING per3"); + AliGenHijing *gener = HijingStandard(); +// impact parameter range + gener->SetImpactParameterRange(11.2, 13.2); + gGener=gener; + } + break; + case kHijing_per4: + { + comment = comment.Append("HIJING per4"); + AliGenHijing *gener = HijingStandard(); +// impact parameter range + gener->SetImpactParameterRange(13.2, 15.); + gGener=gener; + } + break; + case kHijing_per5: + { + comment = comment.Append("HIJING per5"); + AliGenHijing *gener = HijingStandard(); +// impact parameter range + gener->SetImpactParameterRange(15., 100.); + gGener=gener; + } + break; + case kCocktailTRD: + { + comment = comment.Append(" Cocktail for TRD at 5.5 TeV"); + AliGenCocktail *gener = new AliGenCocktail(); + + AliGenParam *phi = new AliGenParam(10, + new AliGenMUONlib(), + AliGenMUONlib::kPhi, + "Vogt PbPb"); + + phi->SetPtRange(0, 100); + phi->SetYRange(-1., +1.); + phi->SetForceDecay(kDiElectron); + + AliGenParam *omega = new AliGenParam(10, + new AliGenMUONlib(), + AliGenMUONlib::kOmega, + "Vogt PbPb"); + + omega->SetPtRange(0, 100); + omega->SetYRange(-1., +1.); + omega->SetForceDecay(kDiElectron); + + AliGenParam *jpsi = new AliGenParam(10, + new AliGenMUONlib(), + AliGenMUONlib::kJpsiFamily, + "Vogt PbPb"); + + jpsi->SetPtRange(0, 100); + jpsi->SetYRange(-1., +1.); + jpsi->SetForceDecay(kDiElectron); + + AliGenParam *ups = new AliGenParam(10, + new AliGenMUONlib(), + AliGenMUONlib::kUpsilonFamily, + "Vogt PbPb"); + ups->SetPtRange(0, 100); + ups->SetYRange(-1., +1.); + ups->SetForceDecay(kDiElectron); + + AliGenParam *charm = new AliGenParam(10, + new AliGenMUONlib(), + AliGenMUONlib::kCharm, + "central"); + charm->SetPtRange(0, 100); + charm->SetYRange(-1.5, +1.5); + charm->SetForceDecay(kSemiElectronic); + + + AliGenParam *beauty = new AliGenParam(10, + new AliGenMUONlib(), + AliGenMUONlib::kBeauty, + "central"); + beauty->SetPtRange(0, 100); + beauty->SetYRange(-1.5, +1.5); + beauty->SetForceDecay(kSemiElectronic); + + AliGenParam *beautyJ = new AliGenParam(10, + new AliGenMUONlib(), + AliGenMUONlib::kBeauty, + "central"); + beautyJ->SetPtRange(0, 100); + beautyJ->SetYRange(-1.5, +1.5); + beautyJ->SetForceDecay(kBJpsiDiElectron); + + gener->AddGenerator(phi,"Phi",1); + gener->AddGenerator(omega,"Omega",1); + gener->AddGenerator(jpsi,"J/psi",1); + gener->AddGenerator(ups,"Upsilon",1); + gener->AddGenerator(charm,"Charm",1); + gener->AddGenerator(beauty,"Beauty",1); + gener->AddGenerator(beautyJ,"J/Psi from Beauty",1); + gGener=gener; + } + break; + case kpieTRD: + { + comment = comment.Append("e, pi for TRD"); + AliGenCocktail *gener = new AliGenCocktail(); + Double_t momen=2.0; + Int_t Npart=200; + AliGenBox *electron = new AliGenBox(Npart); + electron->SetMomentumRange(momen,momen); + electron->SetPhiRange(0,360); + Float_t thmin = EtaToTheta(.9); // theta min. <---> eta max + Float_t thmax = EtaToTheta(-.9); // theta max. <---> eta min + electron->SetThetaRange(thmin,thmax); + electron->SetOrigin(0,0,0); //vertex position + electron->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position + electron->SetPart(11); //GEANT particle type + + AliGenBox *pion = new AliGenBox(Npart); + pion->SetMomentumRange(momen,momen); + pion->SetPhiRange(0,360); + pion->SetThetaRange(thmin,thmax); + pion->SetOrigin(0,0,0); //vertex position + pion->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position + pion->SetPart(211); //GEANT particle type + + AliGenBox *muon = new AliGenBox(Npart); + muon->SetMomentumRange(momen,momen); + muon->SetPhiRange(0,360); + muon->SetThetaRange(thmin,thmax); + muon->SetOrigin(0,0,0); //vertex position + muon->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position + muon->SetPart(13); //GEANT particle type + + AliGenBox *kaon = new AliGenBox(Npart); + kaon->SetMomentumRange(momen,momen); + kaon->SetPhiRange(0,360); + kaon->SetThetaRange(thmin,thmax); + kaon->SetOrigin(0,0,0); //vertex position + kaon->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position + kaon->SetPart(321); //GEANT particle type + + AliGenBox *proton = new AliGenBox(Npart); + proton->SetMomentumRange(momen,momen); + proton->SetPhiRange(0,360); + proton->SetThetaRange(thmin,thmax); + proton->SetOrigin(0,0,0); //vertex position + proton->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position + proton->SetPart(2212); //GEANT particle type + + gener->AddGenerator(electron,"electron",1); + gener->AddGenerator(pion,"pion",1); + gener->AddGenerator(muon,"muon",1); + gener->AddGenerator(kaon,"kaon",1); + gener->AddGenerator(proton,"proton",1); + gGener=gener; + } + break; + default: break; + } + return gGener; +} + +void ProcessEnvironmentVars() +{ + // Run type + if (gSystem->Getenv("CONFIG_RUN_TYPE")) { + for (Int_t iRun = 0; iRun < kRunMax; iRun++) { + if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) { + srun = (PprRun_t)iRun; + cout<<"Run type set to "<Getenv("CONFIG_SEED")) { + sseed = atoi(gSystem->Getenv("CONFIG_SEED")); + } +} diff --git a/TRD/TRDdEdxHistogramsV1.root b/TRD/TRDdEdxHistogramsV1.root new file mode 100644 index 0000000000000000000000000000000000000000..7a96c2fd6ae2a591f46aa17b38ece1d69a3a0ca5 GIT binary patch literal 53434 zcmd43Wl&sex2TJ|OMnCq4#6E7*Wm8%?k>UI-Q6`naCZpq?(Xh($oke<``gLhweLOm z{y0>1RSP|t)#K^+9nTnJVryk(4+PX(1Ox=64+J#c^nP=Ezc2891Af2BgZ=By_fa4q zIMVm0IAt>qYEkQ}iS7?Zz1z;-U;Tgd69okN{(Fj4d=D2u?0+5k{jUK5fy)b-m>Xcq z%L?cV>N|^=*x6ed+3Hx>Dbi8f{yMS${uLmwU&s3EctEJ$?}w6qzx{PgAfT3;|LH_Q zjDMXExXNGeXa4Jau7ABBOkRXe2-=wL*V#ym(!uEq(&#(m;|si>7&X5f6gdcx3OG=? zEQc|g3>v@;MFWGkR}u|E@!4y?0Cjo<=gvasGw%;+Fed|lanytrbkY#o`1zp4Ap8NZ z_;}^}V*9EXZR8ljgwXJi@b$u&^+u~TF+*Lb533u)fV<_ZQyfDfTE`FIodNVh;AFWS zp>!ZU@AgB&*fyR`Rj3tVf;jY@-uqvYHXe21c!!XvUf~%B@ zg+1Arp{X|0D;AKNEuOnS8j+8!9H|9cC2JT^)srJ!S3o&y011A5=EjU$$2fLX5FSI zk!Bz(c`v%?VpPR5>;)EWgip`jFhVSi2OOf38Uj2cZl9@)hKSFG zLRQ<^jQ|;yPOEO!DlakU^vh>MA_*0>dbLKas&9eZKtLdte{rV#ADn@z|IV49_#d3n zfB-4~;tXIcCxZzv`6S^Rl*5nQ8F-R~-#2Rd*t%_agdbC-D@S-8j3%0s7jlmll3FVn zGL%bg2}6y5fJllQKbft;FgqeBeBsV=myoEM&2Pnemv**%wr8sa_P}oIlRbrrtc}g~ z+7RH~wg0SNGkZK+K#Is@mT$-cB|d$cW~AyZ%AfMpg9kR?n>6|iuO&r51q}#Qc1Nqe zhpKsqd;LWCBB@39(cO{Kh%ZC}7*Hd3-FNqk8^epsd&7&# zN~xhmbtb8w$G6Gf#a2V)-cY10mqODp8wPSTx3w*@r$Q`$B-hL(-+I{poR(T~(gG21 z+^zCXzOyHBZm8_OuU4z}99b2ZHPS`T=x*!}%;Y{Mc=kZm>(4AJpk!bktAVdh1rmOE zeAS8qYzpJRL4Lm5ATUS1BqLH{Ew2WaKSTLMH;k!WFjx;}SPXd;bMcdE*cRl9(cfj- zcVmwBYTBVAu2@PYvX3@M=B1AdZ}4oea^MU78}tl{jo*OXs?sLMg|{-XCQ=QUkrK!z z_F0`^nM_UJ(K-ALgO7vGS0ZXpVM7m0)Zs*2l_1#Z;whyxDH^yUUdP518ah4KpO2C- zqUx4kC7kipQ%E4hbhpVc zGE)0_E7LvaI_g^De{Pg9v)0&ym6Cg2Ra5VF{Nb(AYIe+YO^>FfMpF42`7hdH+Rfx-(E8iQeth)xosVEp*F}kw@LvnNu z<`$vD=Y_)gO2&rY{wCCjUe~mT@e@*Y!@N|D3POgIYp|L#)K>MYTP5>0*Cd58m?p9J zPc>j%{4u-P9;{7iamUn^BD)^A7T0D^x(2(Jg|-j9Cx&kM4K0|x>8R$(ePHVz4B=OW zndRuC39~RpT$Nmt8>S2lVJw?cX4evP^o zF};bYRk#}zlO?hf5b1slzH#HQ!I~_f0gXYqK~&|))*!m~mWD+vpGuld;rhhT$V*c3 zFdg$DK+-mjOCk0dM1Q;=7ka#+TmS-V@dtf>h&TiT=AYj@t2%^kbT(G(*Hu}UP zELd>rqZ1p;(ukevk-Ylq0!b2^=jdCp))LWFX8Xn41?>Fl~o)tuxQlvB`~nc>$@DS^Z7^CPVL8p z)ppzWpIfnz8q%9#IV#Oam7gnz?)AArnHGS$sq-k*d zNj-oEalAPt)$ZdSOjC<#88HtpwVul60I7^gnuGgc^?LmbAK}xp9^EtObs^9e&T06Bu%$g7 zT|pmhf~6~ML*Ty~cdW|C=tZ`Y>h?`iM0$^gEq~hYp^k`2nYV2_6tTU&0l=mZGgmu? z8dS{Dr-vp>)kfrsS1KF&Ju^D8li<9`>M4&G0v6*R$@glhFm{VcVwocuNfYv!>|`9b zuPzP95p#~8<+$YZTh&Z%xZir~(l9`|-!n~+)ZdAM3HY5TY1u!BQvVlGvA=dow0GOa!z}uD0SnJse#%B#Z6ozWkRQr$7aXt5R>;jvB z-cG&(>7iYSIW+ng|CBxO{I3pp1~wszX5^~n_U;12jD%GiP5aI|0+wIX8rS-uC*+$< zA5@F07nnKi9H=nV*Wzh1DHVnuF7uoWTI9=*4F??$(1JbcVFMAys7WFUKA;M>akb}< zD8yprEMd><;9?bevgnf>^D@`r=*uEbrejqmi0x5MKc8gWzU)(X53#?wU`EHR*l8?f z$Tby4Jy)h)7$y#PfxDs9NSBocm3@17tsPVQ7^CTSJQAQp!_8!wS?3U4D|w_Z;Og=; zp){VdDZ8}t?QXfhs(L(2T%$*-3&@Y~;D#0+Pe4f48rtxP{u1sYZSL?aj2 zWv542PJgzdEjnHS&J!XczB1?6vr({Vl9SZHRR9&2QZyq3wNK5Rdm)-^C-| zKqARrLC>`z4znBxpp0MqY9LD2Qv;?{%u%KhfNplvr+Q9sjMbUU;EtHA`1Z`i2i%nX zW|vwHzhh0=8LOW)e5v&8My{JS7Pkwr2~w-%vRM+VWeuNE{e+lWS8GVTcQ`$#bG$s+ z-OE_sadJTbq#a znL8Mezq^GQFwiQE*N&XKc%RJXA!kiX#MYN*7xsG(zMY6zeVQ`E38 z6K`B#snM5+Vkk>a?q~Y89jzkkghb&wPC{hu1Ex26R?~}vWCuTIyW20=7;A8M!jnjv zfCSEG_1f6=xI^36$SQ8hxt3;XnQkYGzW}nrQe28B4ubP7vN0GCag8fCqgZK{jPUkQ zTJSAZwd6i;>y89?Xf7>%oO3p1{amZpxeD>}wqtFZGGQ}43^#6i|L13f(L0BSU{5++ zqNv>yNGn5ku%N6sT>J35Gdc5HTZXb*Nk`SElNcXV_n>g=Z~BgZ5k#U*2`n~=fk=TiVJz&x11GzO!) z!i=)z_YxNi!T={fxbTc%^ORHJgL;`WN6`(_cvDeFL4oFafPzAivC`7&(^~_H>U?LF zhB8J6H51eoG^7KSRt7f)1Kr~Oo|k^Jd8A(VbFR_64fhZd3rH8#;Mj|O+Y1v8n}-v- z1Gw2BK2H?_0VN9O7SDPf49&5_QUzUt3KU`qHgM>zqk^{(N9ESkiwGO1wvG|#02Eb^ zN1YXZ7atV?4$fyi2^7Xh7ZjYV!x`cIaRnq3rur!}ZG%3xwl%c1wsjSRE2~Qc6zL+;{)R)DM^-Wj?EN=Ms5BqOJ`Zp6t^868{0ub&^KU!uCGpoa@x1bRtw_{nAXrB$D*S@hK&$WgOXBjd7yR@1hzvGXSP z5twx5N(?3G%yuKsfUiGJ+14ssRJrb|9y!a-+xgW`-56KMN9`&M$XLWPJwd z+dQZW=8NpikS!+9_&xn(BUg~h<&U4U5T&b`bgl^+(<4%{?9U{(m0CH(HECWMVy__-TT0@ z&9|L{V&Ub-JK5#OJ!goN?CN@oBy1MG%CStrIovpuPA9yQwK)1Rh?8n*l?p9t3b)t} zQf}N$VyB4}=@LL=r)=YOd!Q?*RIhU*&NUp>ji*O4`d5H{y!+?n#GNiB$cL3!rInj= z$p%;cp`ccYm{1^6=cfy$D5{J6y25k@4`CFb2ds{C{59_`z~x9F8zn8Lx687h3cR-F z1j~S0M+-qEXc7j(FpNryx5~LrMRZz67!oAy0g^vnTKz`uccG7AnNu{;Dt$@=v{0*H_0&`lOd<`}`)1#*;IcRjD) zY7-k74|0`TMa1MZ+);$MjBixnrxA3{6l%S`H>?%fU~WvYK;L_qidXhjJxb2I9Lt5e z`B&wTI0BN;^DUZ8AInr*gelF+W5hAdGLEM>O>88{D88kpPcPJ2RpG*>z@pv+n#GmB zh@HtxnSMNaTKJVbxc|-_`mb8$f2uc1DEz@4;9uOy0gO>~FgjlOF>1R_aE^6N z$ytq;?uOIX@zILc@fjzea7Cht_)4Gx$(n*#p(21Wc(WC-vBxAgOU4CFF|sXv$b6*T zSI^oi8}H@Vui?k;IM3b^hXltt{L{3DfUluwgRr@Oa7xd6K8SoNtCh^eIVkndR~@?Y zDVpFMt_BIcRQ=kj@aa;kDC8rrC1JoK2k>cDM~(g_oDG?)+Kl2fCM{^8m3TU|C2q+q zjlhU&E!UAoYH&(8K^?XLVw4(F_*Q=g)rWovjL`me8qV>Jl4PvCAN788$*85TB!(3! zp~ojp17(`}YOGUb!igWt7WO5t-0g>YlmKrQzvR54B7M$*G--L-z`dHXni1njLH~fQ z-)X)w9Kdnog8N#ZJ<;sO^&*0sF~^AgiPR11d#u!2EBFH1w`t4zt%ul-VDPMG2q?>X zWcM&+S0C#e65EzS*yU%!oD7l( zlolx$NGErEO;nm-wH~(35dSDARA4%ubJtfX-QGi&!@6()$<@il8L|}g2wOQ@I6SJ& z&uVjz@(Jf&`m2~cRi2Ez_Ow1Yxu+Mu-4@2ewjXj@TZ=Sw1Yl)~+*PV7z(ij9m!0?9e-4w9N_cZWx0z5&WTs4I_u9NBp? z7p7SBUAvXMB#+(ZqpvV6Bpo?b-OP`?Xoc_I;M~7&u^{Tdvxl+%J9}chf3S!1x2w|` zqw%77n4#F=eh5+I)fo&X3!>94B>8rBOTp1-?~{=hRU35izP}Pw3sRNk=<0mNRyHp5 z&)P91j_IRFC?;}b_KClLJ%$^uK%q+88F_kZeS3pzrcBz$M>VNZ;W6DXhMA5MsH;1V zL*o0{aD8TGbjSHbp#B>`^l9bwtFsZgC?x5EE!-B*`4Ho1`nNWP_WAh|t@{fd{y;Nv zqlhUk6vIKvs)YykXXZd|czQFFm=JxsbK zNps62+$0%Dg27F1h^KnFT_g8cFITg3 zz3WX6x_rZ;R+R)r#{S$9%gZmeB7SGL+k z5{wA)kjr}RQqxM41 z>~x7Vp0X~M&PFhN8OlqK`B_WW)MgWpib$`=rSC4?fuqkBiOH>D9+Kf&e5!z$8j$!Q zMDB6cy4be2RN-s9i;?4EsX}D6KGX`&ZX!80;L0B0WN8g>x)t7PV@zcFoO^tmk(8}T z*8C`knM?=Cmc%~Lf}K}St5NE#R66`H2N-!7pw>&Jx8g=yPX^rO-L{as`E`N~ zmwTGO((fsOx0P8mRuoy}6f7(m%fAh@Al$*wLS!ZC+d@!37}Hw_BWRl}9z#Gcev2BI zhX3*`#XU%}(e^+Uua&z*lueukLpi&vNo}Sf6JgtmQ5J{LWqn1k>SN}nYht*J z0%kvPU^c6kJ0DWNDgb-Yl(PT+Wu6&;yhr&6hpgc4Ju}``SgUqc%D(yCh3iF>YMnjg zL=k6`J>)?W5&af_mh>b4XsR~_didy~%{1S>E-OeRU zU>`U@z!(o18cpY1&(C|mN5^&SG%hYxiOI~#VMSh(l6U~8Gonbz&&)|SN7};eOCdSM zO@TFvU||}XV>#7Fis;I^bRvb^X*d&u*<4lb=V6Kh4%%vseIyiC!8o&@@1l-HJ6Xl? zGZZn!mSm@Hr6ZB$12S-!)0?K@RSDjb=6-+Os^wU1T}5MT_B{3!V(|hZfp-9R|M57s zuxkWe>8vS5|5kFl>tkhcd3B}h3l~wA``7`&nBEDr{TDGwb+yFc+3O5pC}r;%jo}7ap*7-_EyY%jREJ8-xP z9}G-&5=PpFA8fP@0-y)`?>Y|5bBwSJ3Hor#a}2!w1oGg`vl!~#*^55lXaS^)zlJ*C z1UV@IPn&x&`TkM=Y3cYx|3u&)OK13REIo$%S9P|p{W`L5UwT3?`<@o(Fh@by!cFsP zb^U`tCvR{xR9k-Buk=+4QmHsL77Wb+X&Ap@kY744Rt;8T)iC`yCsg)jbfnca7oFX? z>ny{r?P`x}HB|eW6UHZ-UI^^HTaIUeNp{-g*n-T$KEptI$GE8%E{P-{TIvyJA;=tV zOu2=_PWv1z^WH8ulq>oIYmFh8AtzMWlaF~a%^_;GOjWkR>BQWGEelow>4pub06G$y9-(g>iXVn!v9X+wWNgVa)He35x$g8_(Za1-CJ(40eC}9QQ~5zz`ooghJ>Xe8mHz z#oNRqSXNfF1Lg10WVpSuF=&+`)?(REgb3g^QCOO00 zL#;`9Us&(E_S)9&o~~Ahno;2-yHT%kdcAfE4+-B8f zdV)4YPjxsY`>3%v{^mccQ|{{X#RU#XbBZ}4)2%uKXQ`k?~GPKmlFs;;s$Z>#DbVX6euoHJ% zdX6{Y8oxR8s`jgM?a!2ND^shHr}P(cnM;Vke9e1p6#%P)4DHKy%f8f=xNyWQUihr?&C2Hs6e;ucjp~ z!>EvEF%HQO@lZ0wFAh;Do@Ug*XR=V?B>D{2Qx%@+-W4kjVm>nm$&`@Y=OL$NEQr)^7vvEX2{ z<4SglTJO8+fd|4JlD%f$osK|4$ef#eyL-+YyW6q#nGn1k{gB zuecsv&e&}^I5n+O#A8<;WE!ePNT^ukua<)@(8#-t;L9Oo-O6|c4qUKVT#DO}3#rv6 z%u%KIvz*ds+Tf-YF(m-nWz5XuPAgnIuZV;rbif9l8GY?AG%faMBrW!g+P}?D@07f^ zH|v<`NUVqz+Bp>8r%EE=)b;sk=yxj~;H~lpGui=DEqH+Z<%vCO)D8xtQ#{tIj11@x zZ$~=rH5@413@p*k46-^NJBk^@h@>r`@MK%;(H zk4!ucp)Koe9Qwg=gS_;EcDk1_oTR`Pa*7B3!!aOiK=Swg!*~Ue1D~q-buS~#qiJtc zxe zm2xUltle|^PVIN*Nd13jPK@^t z=CJ>k^!N6s4jSMfGF(Q$1~o7ceWF(p`((^x7aR;@0uIJUfbjZ#a)UhY>jt_F+7TfR zsGt%gByyrQZ*Msn+Oszt0y+kyMz39)d!@$&8n$b<%Uy=sSnJ)^kP#;JXCqA300K~< z=&eX*sO`JawBm%4ko%2t*KJ)_;q*1L5>$sQC z5sjR3*dZH+`B-z2P2Mc0aWulAq?EO(Ho1z$U%;+vY;rwAzVy?9ZF@HN72B15@*kWC z@9L>kLS$5@qm)a#O|UkU8WDbDSKM9Zjkw0v=Y_-~A7M7VC}(GSYo1z(tP_75uW!df zYEBPyb6lO3{N7MHR<@U25!g6U@d7bA5TN2UR!OSs*T$Z4-n}64NGX;0oIveIo>N3{4i}(0vGN zO&1$!I(btJ>?vK+b*3lWEApWHWa8k_x)sMq;|HY`!8LukYCp+3q4X>aJloRm?S>sQe! z7fYD?%t1~3M2nW+eA+dP0Z8a<_^q{_jip+VGzarZfqQpV(N0 z9dY|unhO*I!_@`LGmgf)rgf{Y7qc+UIXza}ABa7yMa|DmV~_cPJ=~vOEAf{$g-^tt zMp=RV=~Y1co(6KfR}uFc(H=#iUZNjgV7@AQy+WIUw2nS1@q9DoZP$Gcb7k?~O`H(4 zUFE@J#9Z5Makz}&B$^u}Y82tKg(9mP_!JWZei?*m8)=LZ3xj3zfac(}4K<)*o8*_5 z*xnCYJ$w%o&j9N*5Vi$+!+T_tyLgCetldRY_PKi8{9f`3gesA>_P z5>70o&TY~EQ=0gLpb)X7$e}4p{bJRy-P?D^%${4$i~A(pV1tD=rLC4>gRxocxx1RA z#iyXflxX86np*BViaCt<9e`$hxg|Hz9Yf1krNioKrp%`lB{O?XEjjA3TlW#Sygj`M zK9FaK|4PrPO#T7!-!<+%xBLID=K!&8X*7UIX#zq{1Ni_y&SK|RWPxX6EwfE|nUGu& zUx^FRsf0ZL{pn5>8%N}4vv_yn% zL84h;zPQz7NO;_%$McWw0Bu(e${?1J>mKhOZQBP|B>cdgjQ5<|5 z;XBd7k6zdBf{ss}24NQk4G96kYds7M&RgXTf~L(G2nh#g83qqmVf$|AA8@p-lgB&O zV_dp(muW?@DlfJ|1wAzcvSH$xiIdwzWXs6jmb zTH^n$iB8k~PXzwabH?BGoaRe!AI`AEx137agfXN|j)DEpqXX<*GHB#WXHXA52I3 zi0;I#x#-(B)2H_EEVgOTm@o6sDo@Rt7%8W)?^1dU0bN>iX(jkYK_h{Yd6x1lW;Cl9N(c zU0rg0f_@8xK2o?o#t6KXq1xb3yQ1WL)^X4$}3-3RjP+FA!b5O~vYQJ)>51B=k2Amlo8Fs4K6_o!BJ#H*c8tWWk+jmzU5L zW4)s_oV2;JnG7a@wgvvXhA+k8(8zQ7=;L8zqzs#e?3FDHxk#1}JKB{$OtKqYM(=7k zDif{?-KMT;DK+HM>b}iP7p_gcS6xBke$OxHd47jaQ2Y=06#mhrUz2-p2f`@6zPgYX zpt!T0g!3>6K8qQ7ENGp)jMRd((NJ74#3b+IpzH7-mIcs|(OlN0$9cwjZIz zy4ALCnN!+K^LlAS)o!r4X}?|#RlOST^LrOSgVp2sKxm`qUTEWf%*n~V!SsEK!y;Wd z7Gv`lfB}GkiZ)-W176#mUdO5SY2C5*QU;7qY1yv(LNas%3EE#_$(R|Ve)NAMB4I#u7Lg2t~kU*R&zZ`97>;6pN@ zuxRhb1SP)Svz7ZBv1j37=yPO3&XRFZPDzuJaG4%JJ^lVeu&Gt^fZqauJQJO7vD z&(oqKyNo0yQ&-6$b46LiH$9|mHcU4rVur+;W+lBYrA*cd(~M|=hr9)?Kk>!(rWisM z)>+<}^!lAi#_8Xglu-DCN#4J8>b-)bgY99W8lLUjD=2^dT@|fLfsESpO2f9k7GJjB zCs^R(*K!z~vz)`Kl5P(dw2t4v3Ps&1YHv-Ufz3y{qV9)J(RR1uv5^Cg^Xz~3vU>LN z>e0yI8*H!**N3393^X-xAv^@?uQ4Y%)9vbZu=%VHNGq)e&;+HFt0bT+$0?CNATVun z$D-%sD}@la8F_J#ou}2X6C?BrG|xt)KpuC9u_V%f)kn;JOw3ixo!xhyQ_5Am!c;@h z!yb6~5e{pVVg6xE!s^4Ar$ytvsl zCr|f%fO?g~=n}pgrgn(bn#X&yUF*;|N#J-~%SQ{G0@7RCDlF@W0Nc&WI-w2avwEHu zJ(ScGUN)CtR)}^t|9J`^a-i}f7;l(1_a2r)ES8FwA{3}5b0&&4CJ)Eek@E|-o{w4p#JBzm$*t| z^{Sa4!c3QY(2=g<_Q9H)@ii?{k$H-1^^7M`yN1PcyN|=7XUG9Tax@WO_Jg=ECDDUR zo;JDza_N)jmI+59Bjcr-1sm|TBA@0%bxu%hvHMpv$&E$Q; z!LoCZd~D_VWYk#tmddh(roYE}>elf~=l``K|J!vxEPq{}_V3r1iShoy-}{FA_jNzN z<~>FWnBj5IY-9qpIRubPHx(_~#m-1@LhO+RP=nXMWN#A^p>IH>mxrP`uv#O3`_9US z>!Z!Vmds>HK}O97iJW=Y?6yVhXnp=Z?9noI)pmCGi@*y5R)}VVTx~BC)@Qmtug&e7 zSorn{w!MrD)BBom$;M-FohR0%5c8BM#ocBD?PG?T7*^{Kgd*lx>j)w0M&(kWsk$w4 zXJ8uy%^UHlLl`Y;@vcXvD@ry0K2vEZxVxj5kr%|eOp-h5g{Ershh8AJhwB^4fy8^+ zz#Y7u3g8*Uj!i|dN|W{0@=hnhdpYRmK}=5q@4IpctJqrj-6Q0dB@br|;7Z-~DyNIj zdzsOACJu;j4olo+3$dy43goO7}B>_nU5y;PmCFmV1B##2(oCM5J7pth14Czs9 zM?yU?RdENX+qMOB^t&z<_e^}I!x92L0(`9k$Ak8naEGe1iFs{eW4DyT(Do=63M>(|xWq`{J&tDDD@ZT3d zwu#>tKWW)N;FJ7&L;kh+VfJOtE9wo5-$U|Jf~DpVDueU#TC9n+tEp2Fr-GvfM}d+P z2S+LUse!SS$FdrC2*T#%`gU&+NghsiW_J!!=OZ9cZ(eyi&8^m=V?V$5Ywp(W+O$r& zI;Fv}IxF-r(X`R9UK=31j#v8@n$*t3$I^!8RxVr%l^HclJyHx*KAg_=f_V$ z(Z_Pxd!Vi3R=You_(lh>7*$;E!>=bFyV zvXOq%^lg4GkTp+t*cfup=^xDp2TnHo33woKJ!y#A88J-W%~qJF^X$u$objI7Ff41n z4>JA_a7o0Md5AUmv`X$tRX_zu+NSZXZ81h%1x?dt+h}Z(&WS(G@}%0LowmVU1aM;^ zilapXj7#BM*=em0k@R7E%PRjL1=(Jaf14^IyXkj>6LEWTS4rv+VKkn&iIhaMZ?>L4WX|@Xgnz#gB{RQK{_vAS>|AFy8OF!fPx%6BA z-IK3+AD&0@a0>odfkc$KGv%8jxAgL%csfY*ml(4xm=|qq<0H0R(p^>syW#l}OW`K2m3Qx7Er_9LE4gpA z;V~V#_A}Hn@kx235DN9!jYJ0w&Z&P@Qk2h?Z^e$L4!{V6@ic_r^ai#2@v@&8To3et zw+Yy9gkx3r zZiDZG$f#rlDr%m8LJ6Y1RZVVt$^s>FCdM6g^Tv2i^YuV_v@rnGpnLb`XrvQa=*0Ue#ES%5th;1`ZO(^kE(6qj8@73v&-ylQm->wuPq5wBOCxdYywv^ z)KZ~?pv*FGf>|J_IKZwVp*;#@*X!Tj!SVk6^89Pb|F;4Fch5iHvVRHyOurWZxc?~t z{I%pK3F#pfvy%B(#}LZoDfWn3teLfaS}xfyND_Ir!?zHX6p8PG-^|IP-dwj{&xyss z@x}Lp+f=k^ssawA+9W)@@@%Y0TqN+dy7w4+I&)jwihzPqib6dm)8q4^w|$_N;k$8A31LL+N3J;u(o&kif+K1{RqdfQr`;zGRAXD zo=4wtoL@S3#0+POr&C#;vsA1qSnoc(WI2R66|mwpYTfQ0FU%d~U+vsXHlE>S9NOQg zB=0m(SCO0%UA&?zCY;1#e{_jK_lVJ0)g8aKs+(sxm3;ECA6n)aJCOidJ6-0kzG^#t zZ_hhmX`KZ0JQ@H6O9`krSZRe-esNfT9oSQB{hY}@R%wb=yomokTsr2dL2?rsu53lL zh|*xz@vM}!2(JUuagFWlrsphNJcZmc#pI0s2nnXC12k3H>8`(~2DxhkXD6mj>lX((G0MO3wIQ^fn$1U^`ivvamh_oB2!etY5wDjyoDz7`%rzIL&5y9Q>nX@N2s}Gv(T;@Bcy*OFCnXX`gv!?WFnnOq$2=Iu#K(PLtyt|w z-><2DT5WMbzWxx6VsF7838a#DuV)W>5N>oE2fO3CG2ZrSOQ2NYmH}N-t3f4JQzRx< z{AAB}CtMQ7E)pI%?o8PsnD;QOlKCMur(UxJaz!Cge=t3tofN$INtRMxn(Sg_{ zaIX0^3TwjXWYlVFGHRz_x3fcERYPATl|4XyNh<9jMl`H|wnLy5^k&2E3bRhj$IUrK z-EOx|zazM;3@)5T4Ca_l9E=6Cl2|2$2>#MTdv_OnFWK_>Ly*BHNtK8?0WR`YkD+ed zx|s6n*HtW&w@;0knCl$eq|PkQA8b3X9k`fNTK3ZaX<`osf#ccK|+c!D#UMP7X7=EChqJuKLuJc8W~(@1T8m~GTgRKfdUhvPQhbjsWq*Y*jtOkv7Q(ol)t~FG;z`!@#gLpCMmz!CvTmJSAJDTJy7nfm`SJtW0w5XLIg-sG^X_ zGVuGM$CQktEA@ll|4qmsZKK#d4QvWF2?YuNPBd(@^K7O*F><@nYYY`i@E09G*C z52g_3@S*LkvzTI)DieRm@$3Is8DYh7mBc<)Zg8ydcy9aBM8V66R?CrCwOuU~K-fld zn0*qkslkcWmM|@8P%dSgQR#By@to%A$!+!0SI`wroQUsE_A83P5+PO2 zy2Q1-6?kW=@)=<$zVe$lOzSqG@_sf3psvlz5t1b-1d_DH* z3nUV+0yvAYnsxy$^#a~*N0}kbxUUUQ;AE6_>Z|z2S^2|URod_T3W|s=3t&xcQ0HL8 zr>L`17^%8smOV0V)66Y3?B6JzovS9AJ~^m{&sKl?MtAswrb-{q*Q?$N)~KRwkwUA~ zJihRoyL-Wi0t0=wWhQj~Ty+wwe85#?) zwD}zQrk}3wd_8a zP!*{0Xtf@*RwPBO1f<)L+cuk2(ksMVNNHJ*J|aQjK5GS!K}jS+!Q>rz`fYgKYr>5` z>Cg>>*F)*Fpm&oBd({9bGJm zHt8q)fldyrn=S(NaWQpTt<>|z0eVK%a;oqRzPw&sZL0*w*zkpXnbY3jWjfJpgoZkM z!(H;OW`&J<*VQb*?=s{5w9k7Ogu!CA%+UWHx^HGA8(+X{%I;`?9B zy;W4!YrifmB?w5D(hY)icXxM#ba!{Bv~+iabazX4r*tdb`>)e?e)HA$+iS0NvN7mqyfL)g*r_Nx78x~K*<-L(_@3SjL#!w#uX5>)8@_*Uk!pneZ<^Md5tY z1`b?Z*&6eKDq&LtKKE}>E@;VX^)6!6SudjtEY;fuC8cePFiF&qskg6*aa0d2sum{b z#NCfY3oXhi)y4_P6iBfhAHVGR&yqK{GECbICtrwTbREy!ib{`WA5P(<#5T3 z66Sw{@SpER;JWxfX9FyMu>rpxl06u5Py|fSO__p|;FLAr2vZs$(cU4+6$J-1uJcj7 zY6O=mLGg_%%=8#`bVfdXH)IOCVFF7X8`gIWL^!2Y_R zM4Z_W5a~y*uYTFIU~l7?ilBBl8B0phJ19k1KuPN$GV)3UFYd%b_OQeHhShD@s5Mn5 zehORs#ymD}^7P1y`PCkeU8(eWA|vkzAOq-a`?7mX68RmO{_E8Ha-;7?R5@OaLucW< zlZQk!VoRq_^sVoC?B$=}c4EEXTxONO5NuqGmPwH;WW+u2wtaN?;=j!8;nOr9)I`#P zDC5L6%@0#d_~GOhqlfE1y6L4v{^Q!ph%)|{T)=g)|sL4=etipQvBM%uqK<`=WXd)wur9~$Av6N-1Z=sYbX-W5g=QxV0lwi0=>y2T4_xk}jN53w^LIkQ`}0VE{ty4%qJ=s$U@_ zrUeB#mhVY@zR^zjX38ykx0ot0X_*~i;IuJx>1hXga&qMtZf#a_y_Gzn%EbE($bdR` zOaiR*`YHUf_y~8CU{uQxZ}Hc$3DZf>>L8xUD!SH}Otzy+G=H>NXc4l5V=QfOUD1^dagLX)_;74>+RHb!O%4r4-R9;{Zdl zd>N7Ii0R8pn#dMu)vt-12c$_Am4h&>Bj=_2&)s`bhpwzJNvvfQCHt_L6eUS2&}x^l zt%*o_T?JiN@T1u-Xh&EE0m(ClvS{dbf@L!*b@U8WcXgKM?W3P7t7I?O> z@vxiL!4swlXdL$q=F>H1KyO$wP9OM@Z_ZI4RfZ>!Ta7tN-N7X%&gzga|CO%Xr?Q{DQk04pttMW zI1nEZCNUs*_eHT{C%s4>xNRiw5*EOE#XA_tdz+0e&)e)6up*m3&2``CCnzW|GRsc5 zi2L@Sw+>QU;lg!jDJOD1m7nrEpZiY`uM!=Hfij)#19AYjFzQ z;<*7*9$X?){=K3XuCu=B=UJzb{LF-V3HTR!cLf)cxUT6O(3z)`Yjv3gZ?vJaDqgsL z!crXf8+L*oW~Mz=BzruQRl}1^cZpmkm>c{U86)l$|1k=rD#M(xrQ$Kp^R0n5 z4!$Yz&&-o{Gf&hyl5cHq!DXm>oj*pjf;+cA+13-1$5z`yx?lN8En6Hq)MvV9(e1ca zslK*6_tv9g>${fMj~()lYp}w}@9>mV`~#l9uEBs4|K~M$@vR?!C=>y-^b#+&Irx&G z#i;n|5=v-~q%onayB`pecUBAqpSzcIunq8aS$kFuZVi{r&%1X;`D-7kG?EvVE((&PzdG}> zbBL<@n)q?H1vZ<5M1s`~12(LCuw?G+_G)Ecp^t$E=Dls&O|LKT^5uMKH`4Zz2<(CAh zIsv(GOdVLZZkm$r#W}ScnJ%|5`*38Q03_u3gIgn5(6fNrj2m6-fGR)K7`wpX;^3yqL?_7>9w<7j~?WzBji&w z#7&&T&29**&uK*`HN8^}XzkoFMzu9}JfQ6vh|$PSw?I-SH9RtU~18r&_N($<%@+DO=XS_P)6a}o-N zL|cSr=S6LawLGXpEF$Q0=?5_LXFb2>WN3F@h&yGL&83N7UKZ;4;&dxH#mrZAP!iLV zl9VWsTzKi#5V=%Q#U}W_QEx9ZS@Esw9QGH@H(VGNGtC~Bv}^MJ?3!4cFFwq5R-{xo zclYMRSkz=>DNm~sgmooF%xEXyZL;+M_oQBHdn5U>x%3fG!h-tzzDFwC{SHEL-akOd z_tzfzFTLJ;O?YP(fv;Gn4*PU(J`}-B-EQ8g&$?W$udfav;vLzisNASK$W?-dB9=V@ zfzXspzZ`z+cbj-rHc?TXy!F0?DK}~NCuPiPot>DTr>S=Br|E8ur-C#M!5?+z5cMH2 z_yf5kK=A6_T#Attla2h4-QG>;=c$tU+qBwztizJ*HZ|bNg?Y(EKXeF7XZOm!1#KP8 z)?F&uMXHkgWy-dp1}@`kYHe)}^wm_dx$>THwZ=Lv74>k;%Nz~52+FA24``$5i+fu} zeQdW=+S)4avGq0K>9EQG&HoG==W|m*BkqcUrFaHA zJcw@fZd)5?#eRC`jW#pD1jL*K2t*RBoO{p~e&k5nhw`Q7Ul2yeJ+)eE&Sh#ahJMEH z&{}P~hyD?y*DgWfK{^R0rm)NXVk4ZXiwYs!5L0P>XJi|q3}|4q5~nThXN z8=i&WUEirfc9W`CHe=gb`5s>KIXiv%{O0K1>zO=;G@Z1VAdOho}jTD?O|p4I5A+R-R$E zdpx%>(+`JtogqO(yOGE!@IXL38T0d)tbHjm96zF@96hpnDL)did;5W%jaM}!yOqFz z5{F;=O8l{AlOij1O*FLvD8YjM5Ah6P_{R+01$imE0o`$v;m2W~oenYQoF=$74@;1%daeDOc$SnaD=8?O^r#FacG-4#|eW=ze9GYOd955KW= z*R~pm4qJVpbV!xP<8TFbNodzFJhI*_KVvA5dglU(A|>T5GHXI2ft;NKdJ4y!(LD1^ zWED@1uFL6WY4v?`HFYCl#cI*Ly7X$N6~F^5{>}qv*#Er&{yPuA^grhT%zyC!+#;Z| zSWi+T!}!6OS$H_pK-A-=-2T+61q1;i_nL)Tsl8enP4;lts zsJ{`a9VzU7tGJ`W8)yMcE{CpEYBy-v73(NW3>#5rKNab=S`?U=`WgyN>wo!b>Uj}Y zm3z=J;FQ_HGA`)rE_7Wa_?8A$#pfbHNK`b|esMG^L93MDB}T@brsdbBjAgZbQs%W& zs%i7MDba|)_pv6{ijN?047CxxwAiF#{Wbe_947bnMsiZEl7eww@=s3IJf@?!w4?I9 zuh<@x>C}#3a2wPH$5d1rvJ^oB9(4(GF|ObgZ?PR{4_Ih~mvEoxbA(@D#@oZ}oI?)Z z5car0O&s$~55SReXYa+jU~Q^w{iZ#PVa< zGlweSlaWE4e0ho0iQuT5NY!$kAv#L7;@e7uX}h(S ziIq#g04LVd<>yv*DwEXYe=(kaGO7YW+uh;Cj(ziosB~7-&R=n?a}CvSIpYK3=MP zXt&y%l({h9nS6W{xvZ?t>bTx1x1nQGM02;lZE;-AMyO}bEUgk z)`8-&NAU1{^>CVl)n$TFh>?H*CGuSc4(t6( z)vs?;>|H}jRm7L^{m;JUU&1hkLbagxScDEs?zrqO8fMJC^3=I|S-dTn7) zXG zd3jtCv(Vb~x+#abp!&+nLAA76c_&tj)FLgXHBOOEs*F41&eHvdI3Vr!yRNSG??4ok z_ydSq|9;nH2no;P@0B`0FpnnFbwDbuV;+1R%mThowkQJaLBI%jTTM5p5D{|~!P&VA zeSOD^@j`(Z{QWJ0GIa{LV}0~a-rTjuMw3{zFhAypRi~#FyOZWfG66_LIb-(VV9+MS z%}5H^9>>U3GJTktquW^Q0CTIXWD(;E4c5LYjIhPrPghh62wHil2~d1sFROM?W%xgi zX<-@#OCpx-lxW?0TR7kDtV3`{4=0p6$H}!Zq+?)SzT56zp7eu4xU*osx*nnI8MCME zk$%3RAs^{A%+FZiLTa8awVloAl9<5J*B;{S+r~Rym=GJpq^jQ0Hd?%oCI>V7zU)cf zQB8kcH{60x-?`<5-It;%J9O-WCHP)gvO2$S75CsSoxjQMy{Fl?B@yfA?A!^NLk4b? z&q2VWOO8urDB0?pqcw`%mls7Qcpk3@&fE-5^=EiU!S=9{G9SH%ql_-r?#3&S#_jmW zsX-pQc*Z8Zuz0xN_if%U-ctvOUTbbVqcu2a_30lVPZc(|2jK*Ne`ae$6B5mAP$TW)x zub^==3WFAeDO^dRhTvNTZ?U_0iBEskPu?29mLnQgFOI zy+@opVf3t;0paDF@Fyz$Ys+v^%fpZDy>C6HncU>$JJ)N^ne4K|Qz0w!=&)0taOCE%LU!KqB1(T<5TL{llGg1R2se z$Ei+yB}7C>^!<1CMX{P1y+GtN`WihdR?GHv+Yd6JG2vRSbDV+BSMojzSqJ>hlZ}DT zLwrdP7;?>6DtFa;8a>{n^M}t1Y&~@xDOHPKoU1JdJoOZ0Ad^C%nz*2?(OtOuOQpm3 zX{XC;yA3DcTAa=9QoVL{%52EZ8BN1kPw&jd_Zcng&v}^yi1`-fnD@xIJA0O4*A~Ar2mvE#OecNzHOb ztJ8lOr>q<$l`@1}QB^&%C+7}o!N%HJ&S^qoE4VG*e`S0Ac@qASv(N*+>zKyp$3Tc7=JI6Vv%SnTzhd@LO!Gy-CB0k83XXGcBHOB8r zXX=Y8=Mb~D(-AqV=f&_veDNBdkxhQqN~pQRb2zd_!Lzx$CXx0SW>(7;4w0!2wN>&L zt~F~&vh8gC#jG=<(^&1C>W3E`kKEDaB_@=dYsG!EYKKIo$!2AT3U&y^qcwJn(&N`N zTVgv`QWrh9j&(!iv&u8XS4n%}lj#&ow{NeyDW)7b~Q*d{;(&sQq8|y+QfD4eMBCw&D z7J=ph!t_z(1)cpMY-G_pC48tX9j&cn$OtxZ3;ia&!>Hs0c+tPHU#~2+6$S&f4cGUh ztZSUtJFTPxiqVKiY`&#G+*AD+Am81nB**1BrNO5Vo zLi0VkI1^R>3oU~d%;!*!xgW1F{^-MLHrd0^ZyKyyU8x#uN@rUzYEx==LhF}<7t2H~ zlX%bu9%wL)jAogj3;LO%3xL$~oNmQx`XQ4FDCF-VrFUubs&weA3_nBqjOUG+b9&J@ zpI!4^xLn0t#OJt!UPQB=ua76qT}XVs1tbBnk9a9CJ0w-CxFyy{)LL4eG)lhD+f&yG z>MKruMv(fBk$)${_z*SbeF3bb{VX@5RKv*dFlB1A+)}36jmYq0tpS0%LK*1 zuwPUH!6h5~#$7m$@)q-3;36SUtJK-yw+tHUm{TP29)92R2f*N&Ef731c*iX{HcMm( z2XpBvFtHP#XJvP#a0t7PXd~z3F{LNzzZOpy+_udtH^=$bHUEx*S*t2#&~`Fse%5A^ zS*z?CDRKw;b=q#%ErW8!(3CmDeV`sR=BEFU-e{&$p)xSsvoTXWrSwni& zQf&DfdDGjOoHD?zKBHU_cuu-7dTAy4gGIDYbqI*VU&y48yH3EFJawQUV0Z&Qp;4eQ`68x# z@>4{hB=$waESWcn?h@j-9vNntANYPaoEnNCBlH!lReQ;HCix7?Tf}6e(j_<1>O%-{(V47Y0iK{-`*z{6v8ya zek)uiu!>d7VRolTX7#y;%o+1Cm7y@tajwpIa#*l{jC`A`<@1rr+wX=&Rp(dHaC51n zlm!+DBwKMr3qV-ppvIQW$&Ut+))U?F?-li`KWY1%f?TBBP z-U#nS#12-3E*I=h?CSOOR~uR?-flLpEe)_v2HBP$&8$5IthuT!)0{zAoCc6OSxKyK z!l&E(0f|aQ40SB_CIhXQ#EY4N9&z2c24M~&Dny*lW`+I`XjzG?hxlUomaoEe14}(% zJJX8Q%t(~cxX3I1hx(&c=)u$SaIdZ~M+hi*j9br8{VV(0!S_jU9O%V&5hueOE~lC_ zaC%$2+Sqj*p2g%+mY9X?8A&`q3f3tb<)kvylFO-z>2?7k4sk8c!O zEO9UlEPT8B<-iQ;!voD4|H7W%WA;mahw!$HLHYNG1((II17&3N`C}y6C=BAB<4dTU zVEtR+9sqs{Wb=1^!2HwSQ2cj(fcbyU57__W2inEEc~G7vmDK5=Hu1;?eBh4dY$kXu zDlV(YMG!T>F^xk>CDl{VvgHwbd9zr8%{qIz30XpN#5;>DQL!{MW}JLPDBi)foq4*S zZnvG(sFZlTCgI!T?NkkjR%0z;em0GC!41*}*0!pKmMi%Zz zFh#G(5FgH~864RhjSA7BET^a9{k8lUvG&3rzJa@BJ+Nd*F=U1ZW!k2B9IubArg0Ul z^qx-ncS(}0hq70r^(YU+Tzl|X7UAI0S|ySfZV4u4@{&UNaU{pvmuRI>3gA1 z{bCYe5Ry^7b$L9IqfTCa$Mj(kIc9@=sSAABFF+%I?)d*839Rn_UWC6%0+!!N0-3)^ z0y#-wR*T_LS!U@6DF{!1hQ67FIzbk?Db2n#%;QF(Yvn6M&>$lQ6qG$xl72IU5R3$Y zClFlp6Opb84w{UGmV?Kc3=@!U-IM(0QO*8A9+l&U*lO*%^?HOkHjV5;lAz#ok+&8E zO4Qe&2{7kJWjJMtJX%^$U9(=SVg`5ZpxQg$Iv4YUIHjEw1Ku%oLl@a7%n-cFY;X3n z8%i}GAUL~1*3!Qgg#D~5+|KCoz1L@rR+QP2?SoP;2e#6=zy?FyUBc_BaxAV(eMYbR z&&k}6$CsQ*RIF(qBGgID3nAhR+mV5An3xq!(xkXpb6P+Y>)hCv4v1pg@D>d8HPgM| zAp9N2ur+^h!Mx{yC{}q%$!v^}HB>F-x*AL37|ZR(5V12dcFxi>{TX46o0LT4jhID#HVzSJ7?I|b$+kN4=>ed?g$1aCDZxxn*Fgjad$E7 z+*0>@A!Am)WJSm6*%Q`I?;FEPPaeXT_a%`99XUtsYR?7GB~4!qG6znfgdZiPQax^p zvW!$tRn<;mv5!^mdBN;QKaWrYa7_LE0%9BbI~)Zi{(z(AUl)*1Kil#G+b{b2Wb&ST zg84o}9O@q8p-f*9v33ep_P`N9#3J>=#wX{-DEWsaDjv!lsT|2cRt69>X{Pq%S zOqN0U0=~txbQsYRK%_Y4v39!qu$+>DVS*c-7#;t0(omp@xG!kp9nzY|p*^j7sR1>d zJzBL#uNs|+c-#T`;wD3_F@l<{MzxUrBfzt5htmB0bg?&9wETor{yvjDXs+v)et zkEzIL}y$dYV*`cVP2bj~xCKGjqp^P>Vt zzdbR0C<~|8T;mI~>A}gD;iRRk+-vPasPgG6e`%pw%IK!dZ@L0k-!5|MGb@2C0jwuZ zI?9|}LxD+`XbJRX{~4YmT0$oD&{i;qM3^2f=}-@Y9x7z&uzBK5be9fljLI>-O6j@Y z(VgCx1deP7>O8^|SZq2y={)yYnCpVyFc>Qt* z|Ey`R^{RM`e{$Dw{1w6;Xg7VNJ2x!`lI6&U2ex#!qW4{8sbOF(*^Mr11>f$JN{ib$ z=W^cOwG=NX<&^WE!?mOwqS|zq1<*sIKlg7;U*2ofRG5P2kx1{&hS>MJYcp%ruCmni zhcohL3J@NIE1ZY#N{@v62rmMS{vDeRB)?-*Qt=OL%Ke2+`M(;^Ok&>pg@N*EMI>5FXq~2zrF$>lynm~kqsNM)~hDm5~ zBYk(7!V6@qRIm(*^H$r$5cWrd$Z?OXA6=A$KL8e%r=#o0NC&$Sz++`;=KpG6B=wkB z$#jQ+-veEBA`I&+A^CD$$?w-Q?aOjq`F+p}AJ%j0wInCtnf61u9sBq$#cD)S|6H;# z#GpQ&jm^;%mH5bhHQb4ArC?4n!+;Q|ju7MmVgY}c2Og`mUkR@WX4<2YF~6`9S?46Q zx?IOq58{^hKDFxRW`@YH8q3J!=A*h=e^m3oCicJfd88h^hFE3%j?svb-NIjrI2IDb zb4y^8#KH%Qpqkz+h(mNfEGM$3sn$cNo0>>EdLF?&(2dw^@!jkI^p>7gc&BuTdPtL0 z>D6qw((Ds@kDG4#c+-JeX~x-vQxZFPF!=$Fs;*K6!uFz(>BXm!)Y-kmQmBQH`y}y{ zJwfZVLdz>CJ?(+xd1d(fErv_(nu6K9xrqD0#pqRd8~MjJB?Y>e?kRLi_Wz z=2ZVXAjNtA0OYUNTFg&^Kk%8cV#E|tuRvBQ=Z0+9v^Fs@b#PPxtvp7GS6L~>2chc( z8k3PGTxU0DLWdsDS7lGvfK@FdPA9rg7dmV{Pqx(exse{-DBCx{={n{7?v&N(v{Xk2 z4x%%He-B*;3&OME&U^E$t0;f&DBZ5OIeAr}TAixdv)L}Dia?P6T^@HK=3Uz)g~`i+ zq!&$GP}Y!wT%;1|UQWWGkW2(jO|mVHA-KY;~%8%S?aNbxMk+) zFbhaC$`eK-DiZSBTAK?-W{eDyS5}d$x8ibFPXpr&?4nv5{Yv+E16%T>hT2#20ai8n zeuA7nF zft^jpyGPJrwpx3FsUmKhc25+zlM?V=cqhXS-(Vy{?axeh7Gc?y@^M{J&pD2*1+0es zw3}jj!xMU=4Yk3Hc3JVH5ut-qG7zGc`sk94=NPe3oQf8Q5R%m>?CiuyJ7#GMUllJb zH=Rz99h8OJn4)G(RDxiV5Oa(0!^+kzW-7CPd$)^3yl9#Rag!5KxMIQES~dCh?78%H zZ`$bK=SXiR39}k0)8GS<0GCom|jQDx=m zE7v9L)c31e)L9x0I`Y1KE{mL5a-C%ZIYS1Z_(!Xb>mU9<8M!~8_^Vatr=|*nv$8Ba zi%*a@3J(4-lgz^5UEzqwU~>bD$bNoiG?|Auxv0Ma;;P(!KAdkRe@CyRFk~y=e=AAqx>ob_$v{UxAwULpq&T%YgAnZbK{KZ$YGk(k$p>8|Uvlgu$ z9~@nCElCfL@q#xE&3p|6(T^SEixTr$uxk#^!)lh;zTJ8I}Xf4sa_f1dB?JNVoLXbEP_ zj~qy8G4r_KDya1??>YU|!y&a-JzM5my#}{AoAI7$I=EKbykNeQKNUWHpoRv^{Eib0 zttk)*!tM^?s0%%7*_AG?Go(5H6;}tI?m6}b1aD$7}kp}hiq|N0ov^-yW9Y3s5a4! zDcGX%(IFKfI(}a?&efj)2$Gne*X$bVBYnS+8dBg(D;WH&F>}lnx zkc^M;%Qz<3zR5$eap#ri(S0tOwRD7M-iQYcBkIorcu*a+x0YjuevjyCN{__9uh6FB+ z|8th$U)cixvZcQ!_t!y2(dP}zeKjAP31kb*k&7XiQ*sPIpxY%xiJ#T5;S;Ek6UUT@J=jstoQ8!#sD@~-N&9A%un13FdCpo zS(*W#gKWY98rwgs&=xpt?;sTr|$vwN9-@OHo=DJ@QZY!1dq>@Z-8TM1{+ zexyi&<%Ez{Tq|tk307feSn>*XSzpVliqO0!dEzqhhiJA1kJI2Th3mNuJa{;%`Q>O! zj*;(=e8GyVFmuY3l!3mOf{?JMk$&8FW)?TH%F$UN9*4zQ@!{&fa0@ui|~gNc7i(vbNs47BSLW?_I>w5n6@$sG?N~awRzAh?@sl zPMbyPz60ejL0f47ZW86m*~#Y$Nd;)57rViIRf4szK-JCk-yg5Z(#;l?v1>~D%bu5( z6lp)7_w_Fb73A9qii#Q@M$FVw3*>YK>F($c$9>9Na$*zzO3M}caCKF@w&HiZ>OQ5n zg1h2*xo?|1zdTcHB0c8R`c&Ac{`k16S-R)`Hpni^sG*?t@cJ}shH}@bcD~0&)BetS z`g!;VJCOhT3(U9hcMu9n`~gCvzg}Q5f3X8=7s7h>g-S9#_!vRZo$pu_iboQ*z_G#g zj%9i#^99a)KPBfTt3f--eA<-81}_q?d?L`*?CIRt%$zLBoUhdssE330#cNv|Y_-Tn z#NtQ*iUHPG*BW;6#Gfc)3a%8`vh&QXE z8Or*_Pio~Lzil56H3rY}hrzZ3)!<7!bp*5oVcCtT6zsjnj-5}H>6Y}X@f5IZd3@zG zZ?WIliuG!yRv}CB>=mTV&#+b_uwtkOW8GZc1yYqW(Q9fPV^{L0hgdz555q2no0lvM z6k6(A7#1p6OQ8Z&3eMwOvtq;kV zQ_J4b2*Hf$NhI!Jm>aY7ViQl9{3DVMp)~gm3@|^KR=Yk;UyX%|a(|V(xkVAhCBE83 zqYhz+t8>Sh7y^=bD4Y1^fFFG)q=O@mi)3t&vtv1L^_!06M4J#5%rGmup}(u|lZ|b# zCnO^kjXFo-2-PI_0Nt&~D%d{j0mT6rouFH&8i8JJq;;Bo43Wr3bArQH-1L_rXKnJh z1C%cOJFn(MI0=xq-ZhGhMs@4Q=!;U*a)g~uGfncjB{|2Pt5)7J7j%6bK}r}#E_*|w zEE2y=63YdqKoB;~sZbWDC|NhI{E#T2Y}%Xm?f~X|B6Li;-i!v^&1)*IWmW}4yOI%Y z3*BSs@;v=CC;QWkT2>iXCFiuGnML5yXSWMeyG8BYeJBNYc9jaYf=X8h-lJ<)UJr*h ztx_%N5J3MnEEHA9Nn({EL`hTy0|?>z+*5v*<(Eq1lj|bFMz_ev8o61*^+_p!g<_BGp8I} zobA3<=lUDF$1sG0-RJq13yDYk-XK7G0L4K3jjY?T+HN@5b{dZX^|#k&R){dMP%rxs zVM(FBOm7Q<_iKU+?~B7Dm|GF89Ks=a%KBGK>9Ju>_fjn`4lk+{8 zgI@z*tcNQ7HcbPsolpzG4iI3GoFF1rhGUR8leX&JihM~uzjlswF|UY*&^UVC}u@FU(=lbU`WMml48 ziWo)83p84f84neCsatW^dF_qb5A&96g?`Ho&seIuZu@r)yH;Nbop*- zTxOX3yLK!_U5A*gR7>&9BJVC;Ek>@2>RBU_Wh&Ob|46_5!{ZY4qr~Um+9t($|G?nC zJT8B6Rh4Cu{OEDef@Lf!ouT7u2m*JPBuT1`h=X$01c+=YbVP`9f{6jXP}x-WoSANZ z?K;*u4HI!IgX=$U%r$)YfA2u{FIot+& zLq!{hD&;0=EI$?(a>e#8x@};?C(L|6Q0VbeilSxK5O4ubj|tb6dM2EY<$KGVree*z z@H+x$wdRkM3j=OeG?|6S%)OP`LFvmPoaN)aZ_~3uQ~NnX1l_J_DDV(i>r?gZzS%tq zY>6j^p%zqYLdOW}aq`aBAr>9n@txf0o!a*6)GdBqRgljjFT5yyuHOel;wnB45{ONoj;|J?9E@+{m*DBGE(n)$)?OB>>3pFf zTP;(8x!Q;ja57ALU zoB7F@K|OyMt^TkWgew0ILK(R~K={jI@Gqm4Iw1RdOcMC^fxRY^DO=pc;>2tfd)bsc z&7TR$Na>r^ z+aKbpv29}6y(A&24{PN16q&fJo*0@uPwQQc#pO(-1wyVRt__Q@z5@)K)qt#Jky^=i zaPN*!DGJ&)4;%>-y1|blqq}Yt30_dTdXgdh>s7}>%9Q-Xcbz=XnVw5A5AudKp+HL~$Ad^kI(=)X2f zJ?!pkI2ojz`#Pg+fD6Ss*lO(=y6U}WQJSr?4d@5N#2ctj^+*EOy9mc*d-UY?F4%WS znc~yE1Q#Kwd#e@#8T$*ygsvb0PVxeQE}vK=(xqjeI;v-pK3O3r4E6_uoZN_=wMw^; zY9X2R$Zw*_66z*ECQ%~@hj48&oDMXep|X0ZZ?HwPJ^3fKonZsyfk2A0z%XL0K)6hV zA1nU6vc*ov%zpLTX~FDBYmK}&qsce9h8xZpGJ11HhbUhNiM5H5l!uJ3gQ~R_uI&P~ z_7~nLC$QggiMQWv=G3u^9AQbq`hc(36byO3MBhYHx(* z4W~a`zO$p`?q)Kxs4!qwtzFf3EY??^9JZw>#g`9ITarRHOJyMh6a=EZxE}fMve;;B zvB~tX_6VeCbMB|Jv1vBc+b(v51x-lzYtm#Aq$u%UVw~P-I3C+3&Z}UN%xYVWq0|?f z#1L0O*joHl5WH^fU!KQ-G>s&mvS7?Ro?r(;!0U^0ZK>aRT46=0*%?Zl%lsk-?t5^> zU1KX)x!pYTUf!mCmkhEUUsF7zH=}vW5|z6I|C}^|*uggd>CJk?(*lx>%aa=A=C27i zk+}@kgXp&xbIsC2zQOy#hWGLkM}Wuk`9&;OcoJpT5l+x_`wl#?ygyZQ|6~Uu>i*;M z%80W5m+Sz;e~TR;|KG?CRCoO71cXQ__~A1!JT?aA{y^A3gY|zZ1@wVso_YgYzVQ;4 zY7WR6aAcur7?${@A%GmPgxJKvYsc=f0{DQJ-;ih=9IO1IU5c zUxSE|AqUBSBA@U4O<}2i$Ec87 zN-u%P+bNM!{nt>H;)Gn#-?*K2`xX)j1f_0HT$Z0x>4(k1bLOIvVNq!bco12iCIWrj zKV&JuAINT-nwYlWXcP*k8+AGEayYmY3n>gaORxvF)4nr?Xk5TBH3Kch8?*Z|tGcQ1 zQK**8JmWK%)%y```k~c*wL`~wlNq(J)mEW-lZ+4W8IS{ym1*<|T9zo~smo*X@RAqeWmp-5&3J*GHvc)=gx_89t@*9fbC+K7`75gBrCUWa(Y>$77KTxW%Wz1vzI5 z6diz8-Rt@UZphKMOg(xprWIwBE0{khXq$?)vx>7po&69hBkw6dFiE6HCmotQdHsae z3VD+;;0k^JnBl)RJDV|1D&Qf3gEYgjTXHgSm=;A+6q6&sco1SgrE=N=n=I}X*z;^( zytY5_L=^D$@pwstQBByZ&UNF1+f08^{A$ulXU|n{5!XJjYxVJc9_oM8a4L9SP|#<1 zYh-)>F`1M@ZfeDpZ@G2HGQXrqITj^6&9J>Pms6@j>taB~F+YA~^L(6k`hUlvpv0dzbpDH^iv1yGLFw1N%fd%1QxVAR5+f|E ztqZ^buT9I1jrA2=O3^7nd>i$n1cRUyveQqcG&>sywD&*|2z5(W+UMm6fojdiljF@r zuh=}5bUJ(~$Z^wU^?a&lxO#NdBXS;0^pU|{Bob4a8bE+8QA z5ix`^C0Kwb@FRUu%bqzF%!wY{8tR37;-Z@o8S1eY=7AKPCcMj88ToO(F?pWdXWm5n zsO+6xFrX`t!@Gk0g1JC=vkbLl&~qZN6dW6sxvQ-5&ESl)rK>%1G@E{p(H_bD0nr|# zNryW31F|)Gb;kJXJ$h#lT%YS9st&B-4vPjj_)PddpH>6mZfDFU_{B~?hMG>0ad#M_ zOQZ(nB!dR=@>iEWZEr@!?2`KNkCxI-5f;IY`ez1H-ae732|}oKC$`Yuc&g~7dy#IM~E#yq3}AcNl^1Qkj87&Uw2-G#udYc-mx>99g&t6o6=g4f`8YaY=hH9yiAK zE8DqgI%rC2GLyuI^*0434=JWKGf7{+Q5_DKJK@ac*%i)jTfNQzaQjz<_CIhN=k+^o zB^Cd~Es#z4gIxWodI9<~MAb0)f(r32!Cc0nnrdFvP|fZTxmckH9Ffc$gX4<4a`{o8 zflp>BgdnQ{%C?`W*B^v}BAnnzF9V59sR>T{&C7DNHzD5-u*aM{yRX=3X#wOaH9oP3 zw4Q`FsN93Kfw%dIYtaTH+itvA?_ePF(jH61A>{F!Sp(I$3AM%h3xOT>vQ(j{cR}Nq zqTcrdwi%2zf%B+e(_H)EBa#9Wth8I!S2oQzJFHhNcKuI(Tb8H zKizA){xQgfhZ||h`^z4~d@otY_4z%X`r%`)QZjd424xFBBld6_3)?fXtGnpbXVcUz zk`B){54YDbzzD7Tq>z+DOBJm6kmpzss$L-ddfoH&2KH!o%3&9UyqSmdIbutafioN{ z-r1yPLrUqc>NVW$c_kbIu(eNiNeM zryLOkooN0|Y0B1INBK@ss_ep>EeelC#Qm5c0keg$sXNUm(}Mth-bUmz**1DA{2r+RP75>7hL=Ie4SVSK6x1 zXq7Nhlu=``YD$+-LrgWESW7YH4V>G*PV;}DF%hVF|3w}AI?cs-|3o8jn*Ts!?9Y6R z{^y^=+%SGYHV)%X=tZL7n#}y&+EXbu6mDyBVVaklO z>OCD0TUqu2VI!{wVZ+-g%@UB0FDRddHH>eJ8#~uZD42asL|EJXO4oR-xQ1Y#y^8ua9U@Y5$K7TPCGmwN5xHxZ;4yHE`F1$spFBr z57g(TXy*A~>iYxs88DZ0&0`|2Lx+_51(1(iDWJxo+VU9t3{clkN%}mA*Wa|Jgd6am zuh286pnr2hZ#1)!*&S^mBIlM&HBKjviik&#BaO7GdTs1av|MbritHxM*@JlS5>A5* zw+AbcTGz6O1o|h_%e}!3-0VmD9)}B&Iv>-#Wdx&TG0?*Vk}e18w13D0pNwJ`LsGxK=4<-Kd$*IW-)1K(Y#NaXD ze*tCBZ5`xNQNj%uw!_qu8lW5u_zQ()Z{B)hhCUAni0LXsGc%I1jqvc&jf3>`g^~4xZ+#9m9ER=u zY#|pqQyTAar@5X-V(S=dl}{B1AtDdl2xS76R(103S3F2)GiA$3qpO5E>ii-tQQJoM^GWZt0S z2b1CgZy{>)!KR!-gkrOS>f(hvSbd7y=7-l-%KFoeg@m?UFAqGkMr9-6JF$yN9znSJ!SC!P-RFIJi9FeCW!VG;Bs_wsN~sJz3@eH$sI zS9vX>gZZiF+Sd=qh{e46-QbH86RtSH6DbahoeOl8R%F_!x&tmBXmr(2B%^?DMLpP8 zYIobE8&COo9{i9g6N+X++Q^?OKI*Ct9rG1BOGd?m`1sI->%(9g_~RJF8A018)qaAA zw#ZtSZf-SC<^zj(#x(=Ux4cv@v}&urHdN~L7Z?n>_2)DlIMt>pACCwM(3kD~1n2|s z&07*-;r4Q4dULr+CHj?ysOBuP1}{`BY(y<=Z{95L&my>igDpUSl_>oZEF1}*It70Z z2o@4a8fD#-7+5SxPzJA(;c9nNlgGCf$fAj0iI-GiK}`7jkOd31iayN{e+0ENG8I}0 z8!3jRbQ+%emaF0L)e^U-{1CV2XSnD3*aTisLTCgtp-hRpNmZlo^&fQ_v5E3xr=lCM zWU8}tS&xoSi#C5%DQ5kQ{kbo+Sbh2d8={6$RMsI?GTL z#{*`RZD%@25h(0y z)i*d7oqD^3mNDsq93Rf|Psi%L+k3bMxKz$}nzCn~>nD6zIJ2mi2HVa~Da9X*XZE)V zj>h?94{3;t#+gddQBCm40LJjG2vHRoBbZ{IqwKy7%(CC(MhCoRyf?#acoNSxVfBDG zhu+3H`!>$9H*rpf`rA0)++_b>=!9$tu-}agBFU!c5H&S|EjVZ}b)a9jl2hPtJ`>|- zhQ|vW{rc^vFjW{FC-}Y((G{$ufulN|`qWgY9WGWZRMmHeO1pJrHdGhrL=++sQz#yA z)5Y&-8m+`L4^;xpIf@5^kbUZy1$wJ;)$`YL3t@+2G^ggRb|ze6NI3~pTzHb%2oq3~ z1g2?C(-aqDjPM?Si`F&O_*1wZI1O<;F)=hlDD=<1IGmY*>Zpd~bYj*Y zD!SEnEsZGQ{9o+O;3ovrd1z3Sl?)w*CrwxDf`E`8TO5Y=qa$S0ZHN6XI z5bmBM`lH1l+k$gCddgUYb9$M*G0vm9KjnCUUH)lq-&KUiPeorV(_NXqMmn~(#aAdR zSvqRc3y>P(3?(Wn6Sj?ey0f_oA2%90W)2e}D~;i@x@G^gyX!`oV6?uRnwXGQPzFuz$aVKTva(>6$`zzbD@mmr+E2%Z~sCs*0G^s1-NctIGXiR)rg6T9Op%?P=16bgk7Z*8{L&X`D>8; z<5`aQcEi}ic5jxylCiEo?(zq#L0@s~KCh~!KQBBhd^T^z#qYRtT3c7aq{q9hDs?sz z<7hHc6m#MziO?jx3V0jx9$3?PQ`KP!{meR3-~w926?cDxOp!Zch_xWni+e-}=)fk9 z1Xw%M(o^qwCOTeK=QWY-uztuL6hZ~(pLKlpY_MqNMlsH`aT;6W>f5%Wv73m=6e#rAy7da)ps#=PeT5mc?9~q%F4P)PR}6oW&3W27?$8iT^TO- z!n~}?36~jNyds3H6x`DBr10|GJqKbSWg|dSHs??M*Z-WJ0GV(nW&hLB_1{XZ$tq4F ziqO~CK>?=u2sGfZIn`eA^MtWvI7Gl_H)%U~3k>pVxEidvD!7q~cfs$b!sRehvIE8n zY+(u&YEbG!7Zw^kjz@q>En=WjYku$I;_!tl!eh4QSTIjm@_mVj1hw(Wwi&UC>Svk- zl^>Cmef6xqy}aod>c(l{_vHAb{ zxKNMJPgKOd`>NzgPQ~h%vcQscx42biIPMnNrPb3z;**)paK8QRb|VW$!C1)~E zuJi}!vEA1XY%V*NN9n(?Nv?Z1m$zF`Dm_H1`SF?r_cv4^P%>)jgwPk457W_P`sF)h|jrVhOy)r>MRRJrO2gF zp3v1P*}dEII1JmZms^&cd_LPzBLB!K=9*tj2Ky(bX6#PtH8lrV(A5cRC_8I2n_xNQ>*TD(6iwbtos@ z<8X~ep{R!^2w8N_2b}00E(6q-Ds2V|IBdh|J}r^l@|Z}l(N~qW3~6GqPaBuVx{mvh zyr-`=PWbmu)?7pMS1BC)1y=BI^dWe>0hzgs7fnj2+dFLDIZA~=Ywoj${`u3qJ@pVp z2b$7N&-h*mDr0l_&uJr7pFqq>_?b<2-~(RKZ0{ z>xN^NBxawXUirB%`?c^wdf2l6(%OGRk^-ob7Lqtk515VKy@%*6p1?) zIao3*Lf5f+Z(76m1$ry4jzQqs7DCK+4~(64{H1$K}M7 zL|PK+TV(ued8MJ0$%8vN`N<}%&BZvxe+%|;&F4<4^@Pu~aVSzbGK;DTmnZ-d{`kwoPMqu2*3 zq*oFLG>R3QcmkIf<8|ukI?RQvL7tImzt6|#GiOmh%wfkBL{OL^Ljjti4uLA@a$-)Hux%yad~gW1 z;nYqM%c~O{U%pc4Kyd=*{>!7mXvNf*zyT#9F@ zmSkWwXY*ojzXbZa=S=F-8*iNcSk=yi=M4dFtuaP`u;~IrO5M3F-?N0N0iBDwGbyfT zgS|+uvlO=7w32fzAEo-gbgK`B6@Lvwg8i&0;5no}c%8OI5|PH*kT-et*mo$^Epww! z!ktE)e+*C|b~T|5JytW$o0qb(x8Gh`k9v2y-9GANw<8z+C#3NBX*ZeWKTZ@Zim;W`M`C%}`6srUOH zn73aFc8s@U854_?)tH%a$Q+efH=RSuh8GR-dpXXhVwQl;rZ)M{@Dw3AziVbC)!oEr$#5KzQ8+Y-acxi5Go~s+|p1JYW=hd(u$yj1RP5oyMO2*d>&Z!Lqvh_g~+s+!& zW?!6hMYA?_Is0Z4x9Z#sD&3QlxvHD0?%U z^~PJfy;J_i9h9j5I{QGl1D5-92QGi*4p5OT#2$%3YXk#k2P0GYJv6({2>V&+)+i2b zf!Yr14rnDP70eW3h#urZB0j?9X!NsK2$p#|1X}h#7OBv&5jUpX*z9aKkR7|w>w`Mc zH6Oo(Bpli#M}(}XNbM^eo{LcvQ9F*6G1YyV7ilV%@aBOUJ=N|W~r_WfAZVTz*>ma_Z77;r60SkI8EQeaKqKn4Xz%3ts+vO{IcX3?c4 ziU8DsW-0dONEnY8*v^-hc%99oW$6%I&;6bH68EBnhZ>fg(p_Ao?*SQu4?Lqr~y z;fhT&`O>4eG6rSJl56L^K0pyd^ZKWV_2&Ol2Orb-_&Iu?^~vK)gbvzKTfZ&$_)sz; ze-tg!8{)d_-BI zPNhe&hYXfLb)N$6xk2;w$eP!oqyTYn+ZGD(8*xx3@V7w*hy#{?Bo6q%!Ty#wcs-ev zVC{qj6^Io_QAe@tC)xTSg~EXhZ+}*Z%s=&7xSabZdTv;3XjoXu=?}8ZC^(z}7QVxI z<~`yDf)HtSf)Ma*c4m~CFcE)e`MGczI?K}1IjdP$ zS6RG8ZS6z(iu%*_(lBt%yH)pT1RqnDF55Gbv3P z`1E|Vr};}N0r;q6)Ji8B8;XrW{7`B?=d~ryLdbr1i;F!)3170!*J-dznkqC-rgja_f&~P~wPu8ISu-kbm?yA}wSa~Gcn_Wk3bIari-!M7l z12+MIR*|uAbV5>gaxgKlF>#QzHMPSKvvW4Ib|bU4u(dKV2L3oYJywRBM+CE>1{>GZ zvLKT5F;;aapf$u5Oz$Za$cZtX#L<$|J@xebjZ3olaE4to6tUZB| z$;r^c)fgcVIK;5e@ls|8~Kj@9uE1$!1AD!Uf;JL;!)qjhX<13Q4!oYge)zf(BkL>jOgdaKAE!%05JBX#r zu}Hs2*VDpLoHo|~Xadfg0ySfrzMo=+y}f#rmJk!lpQqu9i(XYGxh?28@|sb+RlPW> z;^jUDer+1|r${Qp&NPS2=VomoSWqi*DtL*C)GGAe3KW^(M{4g9{H#>_Tdn-xTM@6> ze(jfkGECQ1xmm~c)(6^(6;|F*(DC{E$bi&@R!iuL&+E~KEk0fvnS#la{7VhqLmSI- z*k?FQg!WzK_Mcx$d$%G1ow98M8(CtLuFHH!HD8>PkWVbaG-a6;GKmkYE0 zaGFpU=v2NlF_6(muM58i&*K=nP=)lY%JUeAS+&&f>-hH?9dNNE;vZ>J27H>?4^ow> z4Y5z;wXa^wL8YFYOG~dQ^%G)DDO>S4SuQYi6aW0EHnVVl&(h&5RP~T$ntJx3oNDxe zw4awAruzApnqBv}_`xsmN$<~`{pVY}i~236s zQV`kIjTl>x08f8ouvM=A`Jy@3`{vYK>vCWCrC6 zlhR1x93BVfLRJj42>dj9)VyCM%V^2_y$QkKBb(73bgxks^5~j|5EG!N*dS&GfJWL#CLLw#Bqy)=gRH z7<7&I{JUkAZr`wdGpJ_k4+h^1D!3e-r5p07<+m6yCU58wP`Aw)N!F^CtWX;g3(-pF zIu&6`?-C1f@5GE_WP>|f~3BC%0sgWPu4RyM5HDN0Qte)q@& z{b{$(dtPWlnfDNvONrOP{%WZG0dG>* z8T4V$N%f_0V(4JCQi?_)x;O7TDomq1PlzUqqS5F8%)+k>+A$2~Ywi!~0^J2!TH$?J zDUCw|rHrew9m@x_!y-fdvk$X|wIa9&h}U1E*w>(RE%nN_QfRlSlA;x$kuMp;7-dcK z&d?eWHp#&Z#fe|*E1w2Ugv=Wa;xJ4ipOi_!Lw6xcbmpB9D;(f*u(CAC?+7`**saHI zefTV){imhlMM+M&L2uXKR_0@rdMTLFdbLOl>zXG0Qr@LS|7DTb;&MFh&yiHl@@R_K zUZZh6^2<3yA2S{z*W)F(reG&jUfU!%LOAXU`MdGp7zH8?t-i)n7EwwUFoALWK&(0Z zm~zF)$SA9s{~KoeqXfTYU)V^}zOiZ?uec0!S0uukG`E3>oZxcBm>(~c-hSAJq4M%l zSdqSR&s$tVUImkVqFNwUHR;_O09pRA?RW|ylJOd1rnY;mYG_R7fWEuifk8!D!3=?Z zEq1AqG_f?PEQaTYdPPbHbqN2AFM?5KX%x7=^h3BH4tqTW`bIU8IAlp%bR3@0 z@2rVg**|rV8tU1gizJEGW<&_t>siz#&BrRVALf1=XD^leHjh@K(`dz_%SuZieNctV z(o?8h8gtT_Rc~95E!}54eVi4%oyZU|#m7KE**T1mzha>vFb}REwoTD#OwQ6VT`aH- zyM|A->do9&le^t$?tm>iePGmJ++6QRn|@R)7gp8Z>a3`88erj)E-bbF%>Tm46Z5Qa z_2p;AsI<&0$hl$X&nEN&;qz~sEnXYtY!M9~Kiu)g%94O7m4^6fj!`@N%A!D7(lB zD0q`3qvp*9Rru++0)tRU-YzICHoo>%W6haK9}f=eEn{tpDtxKdDl!wjDU15Uwje4D z0!yOUd@)zTdLpfG>I3%kp9`c5MGCZFUuxQiv0+<1OQoitu~@o0CsL}i?1jY7dkFfu zid0c2^ulQ^pESru9>LJTk7Wr z+agPA>tN8UDvGxDH4>zWs=m!2qb+#26e+*WnL5j0uLe1Ro;0@8RG{nhpw&QM>S?_y}lEGH{1->`jO`XCoNxI_{@#7+1(kB@hCwAxA zlrCxdDH{jiiQ)$N+qX{$tX-5s(TtxKYf@8o9kYM*^}<-x#}D*P+=1g9*5^&B4|wKbNLb4 zrJ&AC$eLu~ohvxtdT}b#RJZFq*rx1ɳJ_eA43(;saz8&NuH@f#8;AG_(JcCh0 z`I+DQz_`-1>M=nk0XajWbGB3DCE?XZpy}#&UP=u_XS#CPPqU$=@|3J>L6&`uS>^>& z0v<9$wDm1nr5ALK6URPTD*UQp;TA@h3yPE1CBEomDIMV_K7E7wz0(4j{1@u4ofh}yIe?4yty1&v_5LX5e>(_p z(f+^RvwsER4{q=_d}udUN09Kbf;{G3l@p9_peS*fJ19@kCB-%dMain(L1}}ca2G*Q z(js?If0UB`mS6<@peO|;P!#wSX!wZ!Ot8O(kM*836%f8#I`;SQk!1dw+V6uXNV5LV z!$(H->wo>)eQcoN^8sB_iYQQ&xXc}t6DW!r9TX+2dIzNiilPk!MM;a?LERC8rdtF> zDJX%W?uPG?#V^$FD+@f_zZgD7k-r^8L6Yqs!{>FoIKS;aruJX|>sR>LLBkgVx}>ba zpeS*fJE%KNi0osaC|T7zs6S@-+dkq_0YyoR+(Ctb&X8vq6s4d9in^;|!>{`b_1EyR z-#>h}JvzVdKEXc^hTjKKkYxYI@OguVPw3OH3I7V;Uwx#3BR3OtaYXM(L;;?^ph3U0 zq++X}W67%CLEXb&9hl*rrs;q7#|F-||9wN_zx(aGud^=jnD=mq1yJ`e_58af{b#R5 z0Cf*@JOFiPr~k8~8L+bd`#w~_Aa{1(eJzoI$GnHF4}iLd2Md6@hdIWt8UA*7n)!n6 zSx`Fz0Cf-Td0=IiZ%&op-#q_ynQ&l`+k@fv!}GpM)xcEmA-D{n?xBkdpzfhq3ZU); z;-6(D0n|Md0D+aA2Myo9uHOa>@~_E=0gw5|Z^Un}F8{291)%OBTm+!*?7n{%lKC~m z-*(?WOMd{UJBQ~#Yc&8X`@i2|4-E3J^?3u2d3zB54_lQ1*nK)n0oZ*S`+iOGj~u)0 zWBq^Vzye_RX_x_)ck^}&{#K;$?>9IAqx^g2^}v(fr$QWn-KUxtfZeB}5`f*O@C|_7 zUX%ZaRZGD7{(mSn0Y>@vY8`+V05 A2mk;8 literal 0 HcmV?d00001 diff --git a/TRD/histTRDpid.C b/TRD/histTRDpid.C new file mode 100644 index 00000000000..2557a4d7020 --- /dev/null +++ b/TRD/histTRDpid.C @@ -0,0 +1,315 @@ +// //////////////////////////////////////////////////////////////////////////// +// Macro to get TRD signal distribution from all six layers of TRD +// This function takes the data from following directories +// "rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5",// "rmom6", "rmom8", "rmom10" +// corresponding to momenta +// 0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0 +// Running this macro will generate the root file TRDdEdxHistogramsV1.root +// containing histograms for all momenta. +// Prashant Shukla (shukla@physi.uni-heidelberg.de) +// //////////////////////////////////////////////////////////////////////////// + +#if !defined( __CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliTRDtrack.h" +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliESD.h" +#include "AliRun.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#else +#endif + + +//_____________________________________________________________________________ + +Bool_t histTRDpid(Int_t oFile=1) +{ + // + // + + const Int_t kNMom=11; + + TH1F *h1dEdxEL[kNMom]; + TH1F *h1dEdxPI[kNMom]; + TH1F *h1dEdxMU[kNMom]; + TH1F *h1dEdxKA[kNMom]; + TH1F *h1dEdxPR[kNMom]; + TH1F *h1MaxTimBinEL[kNMom]; + TH1F *h1MaxTimBinPI[kNMom]; + Char_t text[200]; + + // Histograms for dedx of e and pi track in TRD + + for (Int_t imom = 0; imom < kNMom; imom++) { + sprintf(text,"h1dEdxEL%01d",imom+1); + h1dEdxEL[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); + h1dEdxEL[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); + h1dEdxEL[imom]->GetYaxis()->SetTitle("Entries"); + h1dEdxEL[imom]->SetLineColor(kRed); + + sprintf(text,"h1dEdxPI%01d",imom+1); + h1dEdxPI[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); + // TH1F *h1dEdxPI[imom] = new TH1F("","",250,0,2500); + h1dEdxPI[imom]->GetXaxis()->SetTitle("dE/dx_{TRD} (arbitrary units)"); + h1dEdxPI[imom]->GetYaxis()->SetTitle("Entries"); + h1dEdxPI[imom]->SetLineColor(kBlue); + + sprintf(text,"h1dEdxMU%01d",imom+1); + h1dEdxMU[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); + h1dEdxMU[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); + h1dEdxMU[imom]->GetYaxis()->SetTitle("Entries"); + h1dEdxMU[imom]->SetLineColor(kGreen); + + sprintf(text,"h1dEdxKA%01d",imom+1); + h1dEdxKA[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); + h1dEdxKA[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); + h1dEdxKA[imom]->GetYaxis()->SetTitle("Entries"); + h1dEdxKA[imom]->SetLineColor(7); + + sprintf(text,"h1dEdxPR%01d",imom+1); + h1dEdxPR[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); + h1dEdxPR[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); + h1dEdxPR[imom]->GetYaxis()->SetTitle("Entries"); + h1dEdxPR[imom]->SetLineColor(6); + + sprintf(text,"h1MaxTimBinEL%01d",imom+1); + h1MaxTimBinEL[imom] = new TH1F(text,"Dist. Time Bin of max. Cluster for e(Red) and pi(Blue)",30,0,30); + h1MaxTimBinEL[imom]->GetXaxis()->SetTitle("Time Bin of Maximum Cluster"); + h1MaxTimBinEL[imom]->GetYaxis()->SetTitle("Entries"); + h1MaxTimBinEL[imom]->SetLineColor(2); + + sprintf(text,"h1MaxTimBinPI%01d",imom+1); + h1MaxTimBinPI[imom] = new TH1F(text,"Time Bin of max. Cluster Pion",30,0,30); + h1MaxTimBinPI[imom]->SetLineColor(4); + } + + // PID + Int_t partCode[5] = + {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + + // Files + + Char_t *FileDir[11] = {"rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5", "rmom6", "rmom8", "rmom10"}; + + Float_t mome[11]={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0}; + + Char_t gAliceFileName[200]; + Char_t esdFileName[200]; + for (Int_t imom = 0; imom < kNMom; imom++) { + + sprintf(gAliceFileName,"~/work/data/%s/galice.root", FileDir[imom]); + sprintf(esdFileName,"~/work/data/%s/AliESDs.root", FileDir[imom]); + + // open run loader and load gAlice, kinematics and header + AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName); + if (!runLoader) { + Error("CheckESD", "getting run loader from file %s failed",gAliceFileName); + return kFALSE; + } + runLoader->LoadgAlice(); + gAlice = runLoader->GetAliRun(); + if (!gAlice) { + Error("CheckESD", "no galice object found"); + return kFALSE; + } + runLoader->LoadKinematics(); + runLoader->LoadHeader(); + + // AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField()); + // open the ESD file + TFile* esdFile = TFile::Open(esdFileName); + if (!esdFile || !esdFile->IsOpen()) { + Error("CheckESD", "opening ESD file %s failed",esdFileName); + return kFALSE; + } + AliESD* esd = new AliESD; + TTree* tree = (TTree*) esdFile->Get("esdTree"); + if (!tree) { + Error("CheckESD", "no ESD tree found"); + return kFALSE; + } + tree->SetBranchAddress("ESD", &esd); + + // loop over events + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + runLoader->GetEvent(iEvent); + AliStack* stack = gAlice->Stack(); + Int_t nParticles = stack->GetNtrack(); + TArrayF vertex(3); + runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex); + Int_t nGene=0, nGenpi=0; + for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) { + TParticle* particle = stack->Particle(iParticle); + if (!particle) continue; + if (particle->Pt() < 0.001) continue; + // if (TMath::Abs(particle->Eta()) > 0.3) continue; + TVector3 dVertex(particle->Vx() - vertex[0], + particle->Vy() - vertex[1], + particle->Vz() - vertex[2]); + if (dVertex.Mag() > 0.0001) continue; + if (TMath::Abs(particle->GetPdgCode()) == partCode[0]) nGene++; + if (TMath::Abs(particle->GetPdgCode()) == partCode[2]) nGenpi++; + } // end loop over particles + printf("Gen Prim. el= %d, Gen Prim. pi= %d in Event No.= %d in File = %s \n", nGene, nGenpi, iEvent, FileDir[imom]); + + // get the event summary data + tree->GetEvent(iEvent); + if (!esd) { + Error("CheckESD", "no ESD object found for event %d", iEvent); + return kFALSE; + } + + /////////////////////////////////////////// + // AliTRDpidESD::MakePID(esd); + // AliESDpid::MakePID(esd); + /////////////////////////////////////////// + + // Initializing number of electrons and pions + Int_t nne=0, nnpi=0; + // loop over tracks + Int_t nTracks = esd->GetNumberOfTracks(); + printf("Found %d tracks. \n",nTracks); + + for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { + AliESDtrack* track = esd->GetTrack(iTrack); + // select tracks of selected particles + Int_t label = TMath::Abs(track->GetLabel()); + if (label > stack->GetNtrack()) continue; // background + TParticle* particle = stack->Particle(label); + // if (TMath::Abs(particle->Eta()) > 0.3) continue; + TVector3 dVertex(particle->Vx() - vertex[0], + particle->Vy() - vertex[1], + particle->Vz() - vertex[2]); + if (dVertex.Mag() > 0.0001) continue; + if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue; + if (track->GetConstrainedChi2() > 1e9) continue; + Double_t p[3]; + track->GetConstrainedPxPyPz(p); + TVector3 pTrack(p); + // Double_t mom = track->GetP(); + + // PID///////////////////////////////////////////////////// + if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; + Int_t iGen = 5; + for (Int_t i = 0; i < 5; i++) { + if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; + } + + // Filling muon kan and proton////////////////////////////// + Int_t Itot=0; + if(track->GetTRDsignal()) Itot=1; + if(!Itot) continue; + + // dE/dx and TEST TRD + Double_t TRDsignal = track->GetTRDsignal(); + // if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal); + // if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal); + // if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal); + // if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal); + // if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal); + + //Filling elcron and pions//////////////////////////////// + const Int_t kNPlane = 6; + // Track Selection + Double_t Ise=1.; + for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) { + Ise*=track->GetTRDsignals(iPlane); + } + if(Ise==0) continue; + + for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) { + Double_t TRDsignal = track->GetTRDsignals(iPlane); + if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal); + if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal); + if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal); + if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal); + if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal); + + if (iGen == 0) h1MaxTimBinEL[imom]->Fill(track->GetTRDTimBin(iPlane)); + if (iGen == 2) h1MaxTimBinPI[imom]->Fill(track->GetTRDTimBin(iPlane)); + } + + // Int_t TRDncls = track->GetTRDncls(); + if (iGen == 0) nne++; + if (iGen == 2) nnpi++; + } // end of loop over tracks + printf("Electron Tracks = %d, Pion Tracks = %d \n", nne, nnpi); + }// end of loop over events + delete esd; + esdFile->Close(); + delete esdFile; + + runLoader->UnloadHeader(); + runLoader->UnloadKinematics(); + delete runLoader; + } // Loop over files + + + Char_t fileprint[200]; + + // draw the histograms if not in batch mode + + new TCanvas; + gPad->SetBorderMode(0); + gPad->SetFillColor(0); + // gStyle->SetStatColor(0); + h1dEdxPI[0]->DrawCopy(); + h1dEdxEL[0]->DrawCopy("SAME"); + // h1dEdxMU[0]->DrawCopy("SAME"); + // h1dEdxKA[0]->DrawCopy("SAME"); + // h1dEdxPR[0]->DrawCopy("SAME"); + + new TCanvas; + gPad->SetBorderMode(0); + gPad->SetFillColor(0); + h1MaxTimBinEL[0]->DrawCopy(); + h1MaxTimBinPI[0]->DrawCopy("SAME"); + + + // write the output histograms to a file + + Char_t fileresponse[200]; + sprintf(fileresponse,"TRDdEdxHistogramsV%d.root",oFile); + + TFile* outputFile = TFile::Open(fileresponse, "recreate"); + if (!outputFile || !outputFile->IsOpen()) { + Error("CheckESD", "opening output file check.root failed"); + return kFALSE; + } + + for (Int_t imom = 0; imom < kNMom; imom++) { + h1dEdxPI[imom]->Write(); + h1dEdxEL[imom]->Write(); + h1dEdxMU[imom]->Write(); + h1dEdxKA[imom]->Write(); + h1dEdxPR[imom]->Write(); + h1MaxTimBinEL[imom]->Write(); + h1MaxTimBinPI[imom]->Write(); + } + + delete [] h1dEdxPI; + delete [] h1dEdxEL; + delete [] h1dEdxMU; + delete [] h1dEdxKA; + delete [] h1dEdxPR; + delete [] h1MaxTimBinEL; + delete [] h1MaxTimBinPI; + delete outputFile; + + // result of check + Info("GetPIDsignals", "GetPIDsignals was successfull"); + return kTRUE; + +} // end of main program + -- 2.43.0