From 17ba220558f67eb0a9a554089bfeda892bfc5356 Mon Sep 17 00:00:00 2001 From: pbatzing Date: Fri, 24 Jan 2014 12:57:34 +0100 Subject: [PATCH] Implementet AliPHOSClusterSelection.cxx and fixed header. --- .../AliPHOSClusterSelection.cxx | 342 +++++++++++++++++- .../AliPHOSClusterSelection.h | 20 +- 2 files changed, 347 insertions(+), 15 deletions(-) diff --git a/PWGGA/PHOSTasks/ClusterSelection/AliPHOSClusterSelection.cxx b/PWGGA/PHOSTasks/ClusterSelection/AliPHOSClusterSelection.cxx index 6c9443b45f7..50447b52320 100644 --- a/PWGGA/PHOSTasks/ClusterSelection/AliPHOSClusterSelection.cxx +++ b/PWGGA/PHOSTasks/ClusterSelection/AliPHOSClusterSelection.cxx @@ -19,6 +19,18 @@ #include "AliInputEventHandler.h" #include "AliMultiInputEventHandler.h" #include "AliVEvent.h" +#include "AliESDEvent.h" +#include "AliAODEvent.h" +#include "AliESDCaloCluster.h" +#include "AliVParticle.h" +#include "AliVTrack.h" +#include "AliAODCluster.h" +#include "AliVCaloCells.h" +#include "AliOADBContainer.h" +#include "AliPHOSGeometry.h" +#include "AliAnalysisManager.h" + + #include "AliPHOSClusterSelection.h" @@ -26,7 +38,7 @@ AliPHOSClusterSelection::AliPHOSClusterSelection() : fMinChargedParticleTrackDistance(-1.), fNotUnfolded(false), fMaxDispR2(-1.), - fMaxDispCoreR2(-1.), + fCoreRadius(-1.), fMaxTOF(-1.) { // Defaults to the most lenient selection allowable @@ -48,9 +60,96 @@ Bool_t AliPHOSClusterSelection::IsSelected(AliVCluster* cluster) const Bool_t IsSelectedCPV(AliVCluster* cluster) const { - if( 0 > SetMinChargedParticleTrackDistance ) + if( 0 > fMinChargedParticleTrackDistance )//No selection on CPV return true; - //TODO: implement positive case + + else{ + Bool_t cpvBit = true; //Changed to false if there is a matching track in the requested radius. + + Double_t dx=cluster->GetTrackDx();//Position for Track matching + Double_t dz=cluster->GetTrackDz(); + + + AliESDEvent* EventESD = dynamic_cast (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD + AliAODEvent* EventAOD = dynamic_cast (AliPHOSClusterSelection::GetCurrentEvent()); + Double_t mf = 0.; // + if(EventAOD) mf = EventAOD->GetMagneticField(); //Positive for ++ and negative for -- + else if(EventESD) mf = EventESD->GetMagneticField(); //Positive for ++ and negative for -- + + + if(EventESD){//if ESD + AliESDCaloCluster * ESDcluster = static_cast (cluster);//Know its ESD so static cast is fine? + TArrayI * itracks = ESDcluster-> GetTracksMatched() ; + if(itracks->GetSize()>0){ + Int_t iTr = itracks->At(0); + if(iTr>=0 && iTrGetNumberOfTracks()){ + AliVParticle* track = fEvent->GetTrack(iTr); + Double_t pt = track->Pt() ; + Short_t charge = track->Charge() ; + Double_t r=AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf) ; + cpvBit=(r>fMinChargedParticleTrackDistance) ; + } + } + } + + + if(EventAOD){//if AOD + AliAODCluster * AODcluster = static_cast (cluster); + int nTracksMatched = AODcluster->GetNTracksMatched(); + if(nTracksMatched > 0) { + AliVTrack* track = dynamic_cast (cluster->GetTrackMatched(0)); + if ( track ) { + Double_t pt = track->Pt(); + Short_t charge = track->Charge(); + Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf) ; + cpvBit=(r>fMinChargedParticleTrackDistance) ; + } + } + } + return cpvBit; + } +} + + +Bool_t IsSelectedUnfolded(AliVCluster* cluster) const +{ + if(!fNotUnfolded) + return true; + else{ + Bool_t NotUnfolded = cluster->GetNExMax()<2;//True if it was not unfolded + return !NotUnfolded; + } +} + +Bool_t IsSelectedDisp(AliVCluster* cluster) const +{ + if(0 > fMaxDispR2) + return true; + else{ + Double_t m02 = 0.,m20 = 0.; + if(fCoreRadius<0){//No core calculation + m02 = cluster->GetM02(); + m20 = cluster->GetM20(); + return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02) ; + } + else{//TODO: Core calculation for general R + AliVCaloCells* cells = static_cast (AliPHOSClusterSelection::GetCurrentEvent()->GetPHOSCells());//Need the cells + AliPHOSClusterSelection::EvalCoreLambdas(cluster, cells, m02, m20); + return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02); + } + } + +} + +Bool_t IsSelectedTOF(AliVCluster* cluster) const +{ + if(0 > fMaxTOF) + return true; + else{ + // Time of Flight (TOF) + Double_t tof = cluster->GetTOF();//Time of Flight for the cluster + return TMath::Abs(tof) < fMaxTOF;//True if the cut is passed + } } AliPHOSClusterSelection* AliPHOSClusterSelection::SetMinChargedParticleTrackDistance(Float_t distance) @@ -61,6 +160,42 @@ AliPHOSClusterSelection* AliPHOSClusterSelection::SetMinChargedParticleTrackDist // being "not applied" or "disabled". fMinChargedParticleTrackDistance = distance; + return this; +} + +AliPHOSClusterSelection* AliPHOSClusterSelection::SetNotUnfolded(Bool_t notUnfolded) +{ + //if notUnfolded true, it rejects Unfolded Clusters + fNotUnfolded = notUnfolded; + return this; +} + +AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxDispR2(Float_t maxR2) +{ + // 'maxR2' sets the maximum allowed dispersion. + // If 'maxR2' is negative, then all clusters are selected, the selection + // being "not applied" or "disabled". + + fMaxDispR2 = maxR2; + return this; +} + +AliPHOSClusterSelection* AliPHOSClusterSelection::SetCoreRadius(Float_t radius) +{ + // 'radius' sets the core radius used for the dispersion cut evaluation. + // If 'radius' is negative, the whole cluster is used. + + fCoreRadius = radius; + return this; +} + +AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxTOF(Float_t maxTOF) +{ + // 'maxTOF' sets the maximum allowed time of flight for the cluster. + // If 'maxTOF' is negative, all clusters are selected and the selection is "disabled". + + fCoreRadius = radius; + return this; } TString AliPHOSClusterSelection::ToString() const @@ -70,25 +205,214 @@ TString AliPHOSClusterSelection::ToString() const // is limited by the precision given in the formatting of the string. // Take care that the precision is sufficient for your needs. - return TString::Format("%f_%i_%f_%f_%f", + return TString::Format("%.1f_%i_%.1f_%.1f_%.1f", fMinChargedParticleTrackDistance, fNotUnfolded, fMaxDispR2, - fMaxDispCoreR2, + fCoreRadius, fMaxTOF ); } -Float_t AliPHOSClusterSelection::SetMinChargedParticleTrackDistance(const TString& string) +Float_t AliPHOSClusterSelection::GetMinChargedParticleTrackDistance(const TString& string) +{ + //TObjArray * objarray = string.Tokenize("_"); + TObjString * objstring = string.Tokenize("_")->At(0); + Float_t flt= objstring->GetString()->Atof(); + delete objstring; + return flt; +} + +Bool_t AliPHOSClusterSelection::GetUnfolded(const TString& string) +{ + TObjString * objstring = string.Tokenize("_")->At(1); + Bool_t bl = objstring->GetString()->Atoi(); + delete objstring; + return bl; +} + +Float_t AliPHOSClusterSelection::GetMaxDispR2(const TString& string) +{ + TObjString * objstring = string.Tokenize("_")->At(2); + Float_t flt = objstring->Atof(); + delete objstring; + return flt; +} + +Float_t AliPHOSClusterSelection::GetCoreRadius(const TString& string) +{ + TObjArray * objstring = string.Tokenize("_")->At(3); + Float_t flt = objstring->Atof(); + delete objstring; + return flt; +} + +Float_t AliPHOSClusterSelection::GetMaxTOF(const TString& string) { - TObjArray * objarray = string.Tokenize("_"); - Float_t flt = objarray->At(0)->Atof(); - delete objarray; + TObjArray * objstring = string.Tokenize("_")->At(4); + Float_t flt = objstring->Atof(); + delete objstring; return flt; } +Double_t AliPHOSClusterSelection::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, Double_t mf){ + //Parameterization of LHC10h period + //_true if neutral_ + //Copied from Pi0Flow task + + Double_t meanX=0; + Double_t meanZ=0.; + Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+ + 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219); + Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ; + + if(mf<0.){ //field -- + meanZ = -0.468318 ; + if(charge>0) + meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ; + else + meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ; + } + else{ //Field ++ + meanZ= -0.468318; + if(charge>0) + meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ; + else + meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ; + } + + Double_t rz=(dz-meanZ)/sz ; + Double_t rx=(dx-meanX)/sx ; + return TMath::Sqrt(rx*rx+rz*rz) ; +} + +//____________________________________________________________________________ +void AliPHOSClusterSelection::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells * cells,Double_t &m02, Double_t &m20){ + //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum + //Copied from pi0flowtask + + const Double_t rCut=fCoreRadius;//TODO: This now works only for R=4.5 as the dispersion cut is not tuned for other values + const Double32_t * elist = clu->GetCellsAmplitudeFraction() ; +// Calculates the center of gravity in the local PHOS-module coordinates + Float_t wtot = 0; + const Int_t mulDigit=clu->GetNCells() ; + Double_t xc[mulDigit] ; + Double_t zc[mulDigit] ; + Double_t wi[mulDigit] ; + Double_t x = 0 ; + Double_t z = 0 ; + const Double_t logWeight=4.5 ; + for(Int_t iDigit=0; iDigitGetCellAbsId(iDigit) ; + + AliOADBContainer geomContainer("phosGeo");//Initialize Geometry + geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes"); + TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes"); + AliPHOSGeometry * fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ; + for(Int_t mod=0; mod<5; mod++) { + if(!matrixes->At(mod)) { + if( fDebug ) + AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo)); + continue; + } + + + fPHOSGeo->AbsToRelNumbering(absId, relid) ; + fPHOSGeo->RelPosInModule(relid, xi, zi); + xc[iDigit]=xi ; + zc[iDigit]=zi ; + Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ; + wi[iDigit]=0. ; + if (clu->E()>0 && ei>0) { + wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ; + Double_t w=wi[iDigit]; + x += xc[iDigit] * w ; + z += zc[iDigit] * w ; + wtot += w ; + } + } + if (wtot>0) { + x /= wtot ; + z /= wtot ; + } + + wtot = 0. ; + Double_t dxx = 0.; + Double_t dzz = 0.; + Double_t dxz = 0.; + Double_t xCut = 0. ; + Double_t zCut = 0. ; + for(Int_t iDigit=0; iDigit0.) { + Double_t xi= xc[iDigit] ; + Double_t zi= zc[iDigit] ; + if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){ + xCut += w * xi ; + zCut += w * zi ; + dxx += w * xi * xi ; + dzz += w * zi * zi ; + dxz += w * xi * zi ; + wtot += w ; + } + } + + } + if (wtot>0) { + xCut/= wtot ; + zCut/= wtot ; + dxx /= wtot ; + dzz /= wtot ; + dxz /= wtot ; + dxx -= xCut * xCut ; + dzz -= zCut * zCut ; + dxz -= xCut * zCut ; + + m02 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; + m20 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; + } + else { + m20=m02=0.; + } + +} + +//_____________________________________________________________________________ +Bool_t AliAnalysisTaskPi0Flow::TestLambda(Double_t pt,Double_t l1,Double_t l2){ + //Evaluates if lambdas correspond to photon cluster + //Tuned using pp data + //copied from Pi0FlowTask + Double_t l2Mean, l1Mean, l2Sigma, l1Sigma, c, R2; + if(fCoreRadius<0){ + l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ; + l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ; + l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt; + l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt; + c=-0.35-0.550*TMath::Exp(-0.390730*pt) ; + R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + + 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma + + 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ; + } + else{ + //For core radius R=4.5 + l1Mean = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ; + l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ; + l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ; + l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ; + c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ; + R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + + 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma + + 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ; + } + return (R2