1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 #include "AliAnalysisManager.h"
19 #include "AliInputEventHandler.h"
20 #include "AliMultiInputEventHandler.h"
21 #include "AliVEvent.h"
22 #include "AliESDEvent.h"
23 #include "AliAODEvent.h"
24 #include "AliESDCaloCluster.h"
25 #include "AliVParticle.h"
26 #include "AliVTrack.h"
27 #include "AliAODCluster.h"
28 #include "AliVCaloCells.h"
29 #include "AliOADBContainer.h"
30 #include "AliPHOSGeometry.h"
31 #include "AliAnalysisManager.h"
35 #include "AliPHOSClusterSelection.h"
37 AliPHOSClusterSelection::AliPHOSClusterSelection()
38 : fMinChargedParticleTrackDistance(-1.),
44 // Defaults to the most lenient selection allowable
48 AliPHOSClusterSelection::~AliPHOSClusterSelection()
52 Bool_t AliPHOSClusterSelection::IsSelected(AliVCluster* cluster) const
54 return IsSelectedCPV(cluster)
55 && IsSelectedUnfolded(cluster)
56 && IsSelectedDisp(cluster)
57 && IsSelectedDispCore(cluster)
58 && IsSelectedTOF(cluster);
61 Bool_t IsSelectedCPV(AliVCluster* cluster) const
63 if( 0 > fMinChargedParticleTrackDistance )//No selection on CPV
67 Bool_t cpvBit = true; //Changed to false if there is a matching track in the requested radius.
69 Double_t dx=cluster->GetTrackDx();//Position for Track matching
70 Double_t dz=cluster->GetTrackDz();
73 AliESDEvent* EventESD = dynamic_cast<AliESDEvent*> (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD
74 AliAODEvent* EventAOD = dynamic_cast<AliAODEvent*> (AliPHOSClusterSelection::GetCurrentEvent());
76 if(EventAOD) mf = EventAOD->GetMagneticField(); //Positive for ++ and negative for --
77 else if(EventESD) mf = EventESD->GetMagneticField(); //Positive for ++ and negative for --
81 AliESDCaloCluster * ESDcluster = static_cast<AliESDCaloCluster*> (cluster);//Know its ESD so static cast is fine?
82 TArrayI * itracks = ESDcluster-> GetTracksMatched() ;
83 if(itracks->GetSize()>0){
84 Int_t iTr = itracks->At(0);
85 if(iTr>=0 && iTr<fEvent->GetNumberOfTracks()){
86 AliVParticle* track = fEvent->GetTrack(iTr);
87 Double_t pt = track->Pt() ;
88 Short_t charge = track->Charge() ;
89 Double_t r=AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf) ;
90 cpvBit=(r>fMinChargedParticleTrackDistance) ;
97 AliAODCluster * AODcluster = static_cast<AliAODCluster*> (cluster);
98 int nTracksMatched = AODcluster->GetNTracksMatched();
99 if(nTracksMatched > 0) {
100 AliVTrack* track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
102 Double_t pt = track->Pt();
103 Short_t charge = track->Charge();
104 Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf) ;
105 cpvBit=(r>fMinChargedParticleTrackDistance) ;
114 Bool_t IsSelectedUnfolded(AliVCluster* cluster) const
119 Bool_t NotUnfolded = cluster->GetNExMax()<2;//True if it was not unfolded
124 Bool_t IsSelectedDisp(AliVCluster* cluster) const
129 Double_t m02 = 0.,m20 = 0.;
130 if(fCoreRadius<0){//No core calculation
131 m02 = cluster->GetM02();
132 m20 = cluster->GetM20();
133 return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02) ;
135 else{//TODO: Core calculation for general R
136 AliVCaloCells* cells = static_cast<AliVCaloCells*> (AliPHOSClusterSelection::GetCurrentEvent()->GetPHOSCells());//Need the cells
137 AliPHOSClusterSelection::EvalCoreLambdas(cluster, cells, m02, m20);
138 return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02);
144 Bool_t IsSelectedTOF(AliVCluster* cluster) const
149 // Time of Flight (TOF)
150 Double_t tof = cluster->GetTOF();//Time of Flight for the cluster
151 return TMath::Abs(tof) < fMaxTOF;//True if the cut is passed
155 AliPHOSClusterSelection* AliPHOSClusterSelection::SetMinChargedParticleTrackDistance(Float_t distance)
157 // 'distance' set the minimal allowable distance between the cluster
158 // and the nearest extrapolated track.
159 // If 'distance' is negative, then all clusters are sellected, the selection
160 // being "not applied" or "disabled".
162 fMinChargedParticleTrackDistance = distance;
166 AliPHOSClusterSelection* AliPHOSClusterSelection::SetNotUnfolded(Bool_t notUnfolded)
168 //if notUnfolded true, it rejects Unfolded Clusters
169 fNotUnfolded = notUnfolded;
173 AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxDispR2(Float_t maxR2)
175 // 'maxR2' sets the maximum allowed dispersion.
176 // If 'maxR2' is negative, then all clusters are selected, the selection
177 // being "not applied" or "disabled".
183 AliPHOSClusterSelection* AliPHOSClusterSelection::SetCoreRadius(Float_t radius)
185 // 'radius' sets the core radius used for the dispersion cut evaluation.
186 // If 'radius' is negative, the whole cluster is used.
188 fCoreRadius = radius;
192 AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxTOF(Float_t maxTOF)
194 // 'maxTOF' sets the maximum allowed time of flight for the cluster.
195 // If 'maxTOF' is negative, all clusters are selected and the selection is "disabled".
197 fCoreRadius = radius;
201 TString AliPHOSClusterSelection::ToString() const
203 // returns a string an quasi-unique string for whatever selection
204 // parameters the instance contains. The uniqueness of the string
205 // is limited by the precision given in the formatting of the string.
206 // Take care that the precision is sufficient for your needs.
208 return TString::Format("%.1f_%i_%.1f_%.1f_%.1f",
209 fMinChargedParticleTrackDistance,
218 Float_t AliPHOSClusterSelection::GetMinChargedParticleTrackDistance(const TString& string)
220 //TObjArray * objarray = string.Tokenize("_");
221 TObjString * objstring = string.Tokenize("_")->At(0);
222 Float_t flt= objstring->GetString()->Atof();
227 Bool_t AliPHOSClusterSelection::GetUnfolded(const TString& string)
229 TObjString * objstring = string.Tokenize("_")->At(1);
230 Bool_t bl = objstring->GetString()->Atoi();
235 Float_t AliPHOSClusterSelection::GetMaxDispR2(const TString& string)
237 TObjString * objstring = string.Tokenize("_")->At(2);
238 Float_t flt = objstring->Atof();
243 Float_t AliPHOSClusterSelection::GetCoreRadius(const TString& string)
245 TObjArray * objstring = string.Tokenize("_")->At(3);
246 Float_t flt = objstring->Atof();
251 Float_t AliPHOSClusterSelection::GetMaxTOF(const TString& string)
253 TObjArray * objstring = string.Tokenize("_")->At(4);
254 Float_t flt = objstring->Atof();
260 Double_t AliPHOSClusterSelection::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, Double_t mf){
261 //Parameterization of LHC10h period
263 //Copied from Pi0Flow task
267 Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
268 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);
269 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) ;
271 if(mf<0.){ //field --
274 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)) ;
276 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)) ;
281 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)) ;
283 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)) ;
286 Double_t rz=(dz-meanZ)/sz ;
287 Double_t rx=(dx-meanX)/sx ;
288 return TMath::Sqrt(rx*rx+rz*rz) ;
291 //____________________________________________________________________________
292 void AliPHOSClusterSelection::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells * cells,Double_t &m02, Double_t &m20){
293 //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
294 //Copied from pi0flowtask
296 const Double_t rCut=fCoreRadius;//TODO: This now works only for R=4.5 as the dispersion cut is not tuned for other values
297 const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
298 // Calculates the center of gravity in the local PHOS-module coordinates
300 const Int_t mulDigit=clu->GetNCells() ;
301 Double_t xc[mulDigit] ;
302 Double_t zc[mulDigit] ;
303 Double_t wi[mulDigit] ;
306 const Double_t logWeight=4.5 ;
307 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
311 Int_t absId = clu->GetCellAbsId(iDigit) ;
313 AliOADBContainer geomContainer("phosGeo");//Initialize Geometry
314 geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
315 TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
316 AliPHOSGeometry * fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
317 for(Int_t mod=0; mod<5; mod++) {
318 if(!matrixes->At(mod)) {
320 AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
325 fPHOSGeo->AbsToRelNumbering(absId, relid) ;
326 fPHOSGeo->RelPosInModule(relid, xi, zi);
329 Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
331 if (clu->E()>0 && ei>0) {
332 wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
333 Double_t w=wi[iDigit];
334 x += xc[iDigit] * w ;
335 z += zc[iDigit] * w ;
350 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
351 Double_t w=wi[iDigit];
353 Double_t xi= xc[iDigit] ;
354 Double_t zi= zc[iDigit] ;
355 if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
376 m02 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
377 m20 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
385 //_____________________________________________________________________________
386 Bool_t AliAnalysisTaskPi0Flow::TestLambda(Double_t pt,Double_t l1,Double_t l2){
387 //Evaluates if lambdas correspond to photon cluster
388 //Tuned using pp data
389 //copied from Pi0FlowTask
390 Double_t l2Mean, l1Mean, l2Sigma, l1Sigma, c, R2;
392 l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
393 l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
394 l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
395 l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
396 c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
397 R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
398 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
399 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
402 //For core radius R=4.5
403 l1Mean = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
404 l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
405 l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
406 l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
407 c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
408 R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
409 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
410 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
412 return (R2<fMaxDispR2*MaxDispR2);//TODO: Check if this should be 2. order
416 AliVEvent* AliPHOSClusterSelection::GetCurrentEvent() const
418 // Hackish way of getting the current event.
419 // Its probably not appropriate to call this function outside of
420 // AliAnalysisTaskSE::UserExec
422 AliAnalysisManager* analysisManager = dynamic_cast<AliAnalysisManager*>(AliAnalysisManager::GetAnalysisManager());
423 AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(analysisManager->GetInputEventHandler());
424 AliMultiInputEventHandler *multiInputHandler = dynamic_cast<AliMultiInputEventHandler *>(fInputHandler);
425 if (multiInputHandler)
426 inputHandler = dynamic_cast<AliInputEventHandler *>(multiInputHandler->GetFirstInputEventHandler());
428 AliVEvent* inputEvent = dynamic_cast<AliVEvent*>(inputHandler->GetEvent());
430 AliError("Was not able to retrieve event!");