PWGCF/FEMTOSCOPY/UNICOR/AliUnicorEventAliceESD.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2008
  17
  18 //=============================================================================
  19 // AliUnicorEvent-AliESD interface                                                   //
  20 //=============================================================================
  21
  22 #include <cmath>
  23 #include "TFile.h"
  24 #include "TH1.h"
  25 #include "TGraph.h"
  26 #include "AliESDVZERO.h"
  27 #include "AliESDEvent.h"
  28 #include "AliMultiplicity.h"
  29 #include "AliUnicorEventAliceESD.h"
  30
  31 ClassImp(AliUnicorEventAliceESD)
  32
  33 //=============================================================================
  34 AliUnicorEventAliceESD::AliUnicorEventAliceESD(AliESDEvent *esd) : AliUnicorEvent(), fViper(0), fESD(esd)//, fPhysicsSelection(0)
  35 {
  36   // constructor
  37
  38   //  printf("%s object created\n",ClassName());
  39   if (!fESD) fESD = new AliESDEvent();
  40   //  fPhysicsSelection = new AliPhysicsSelection();
  41
  42   // V0 percentile graph for centrality determination
  43
  44   //  TFile::Open("$ALICE_ROOT/ANALYSIS/macros/AliCentralityBy1D_137161_v4.root","read");
  45   TFile::Open("$ALICE_ROOT/ANALYSIS/macros/AliCentralityBy1D_137366_v2.root","read");
  46   const TH1F *hi = (const TH1F*) gFile->Get("hmultV0_percentile");
  47   //const TH1D *hi = (const TH1D*) gFile->Get("hNtracks_percentile");
  48   //const TH1F *hi = (const TH1F*) gFile->Get("hNclusters1_percentile");
  49   fViper = new TGraph((const TH1*) hi);
  50   gFile->Close();
  51 }
  52 //=============================================================================
  53 AliUnicorEventAliceESD::~AliUnicorEventAliceESD() {
  54
  55   // destructor
  56
  57 }
  58 //=============================================================================
  59 Bool_t AliUnicorEventAliceESD::Good() const
  60 {
  61   // event cuts
  62
  63   //  if (!fPhysicsSelection->IsCollisionCandidate(fESD)) return kFALSE;
  64   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
  65   if (!vtx->GetStatus()) return kFALSE;
  66   if (fabs(Zver())>1) return kFALSE;
  67   if (NGoodParticles()<9) return kFALSE;
  68   return kTRUE;
  69 }
  70 //=============================================================================
  71 Double_t AliUnicorEventAliceESD::Centrality() const {
  72
  73   // centrality between 0 (central) and 1 (very peripheral)
  74
  75   // V0 multiplicity, not linearized
  76
  77   AliESDVZERO *v0 = fESD->GetVZEROData();
  78   float multa = v0->GetMTotV0A();
  79   float multc = v0->GetMTotV0C();
  80   double cent = fViper->Eval(multa+multc)/100.0;
  81
  82   // number of tracks
  83
  84   // double cent = fViper->Eval(((AliVEvent*) fESD)->GetNumberOfTracks())/100.0;
  85
  86   // number of clusters in second layer of SPD
  87
  88   //  double nspd2 = fESD->GetMultiplicity()->GetNumberOfITSClusters(1);
  89   //  double zv = fESD->GetPrimaryVertexSPD()->GetZ();
  90   //  const double pars[] = {8.10030e-01,-2.80364e-03,-7.19504e-04};
  91   //  zv -= pars[0];
  92   //  float corr = 1 + zv*(pars[1] + zv*pars[2]);
  93   //  nspd2 = corr>0? nspd2/corr : -1;
  94   //  double cent = fViper->Eval(nspd2)/100.0;
  95
  96   return cent;
  97 }
  98 //=============================================================================
  99 Bool_t AliUnicorEventAliceESD::ParticleGood(Int_t i, Int_t pidi) const
 100 {
 101   // track cuts and particle id cuts; pidi=0 means take all species
 102   // consider using the standard ESDcut
 103
 104   // track quality cuts
 105
 106   AliESDtrack *track = fESD->GetTrack(i);
 107   if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0;        // TPC refit
 108   if (!track->IsOn(AliESDtrack::kITSrefit)) return 0;        // ITS refit
 109   if (track->GetTPCNcls() < 90) return 0;                    // number of TPC clusters
 110   if (track->GetKinkIndex(0) > 0) return 0;                  // no kink daughters
 111   const AliExternalTrackParam *tp = GetTrackParam(i);
 112   if (!tp) return 0;                                         // track param
 113   if (fabs(tp->Eta())>0.8) return 0;                         // fiducial pseudorapidity
 114
 115   //  double pi9 = TMath::Pi()/9.0;
 116   //  double eta = tp->Eta();
 117   //  double phi = ParticlePhi(i);
 118   //  if (eta>0 && phi>-8*pi9 && phi<-7*pi9) return 0; // A10
 119   //  if (eta>0 && phi> 1*pi9 && phi< 2*pi9) return 0; // A01
 120   //  if (tp->Pt()<0.2) return 0;                      // lower pt cutoff
 121
 122   Float_t r,z;
 123   track->GetImpactParametersTPC(r,z);
 124   if (fabs(z)>1.0) return 0;                          // impact parameter in z
 125   if (fabs(r)>1.0) return 0;                          // impact parameter in xy
 126   if (r==0) return 0;
 127
 128   //TBits shared = track->GetTPCSharedMap();
 129   //if (shared.CountBits()) return 0;                 // no shared clusters; pragmatic but dangerous
 130
 131   if (pidi==0) return 1;
 132
 133   // pid
 134
 135   if (!track->IsOn(AliESDtrack::kTPCpid)) return 0;
 136   Double_t p[AliPID::kSPECIES];
 137   track->GetTPCpid(p);
 138   Int_t q = tp->Charge();
 139
 140   if (pidi == -211) return p[AliPID::kPion]+p[AliPID::kMuon]>0.5 && q==-1;
 141   else if (pidi == 211) return p[AliPID::kPion]+p[AliPID::kMuon]>0.5 && q==1;
 142   else if (pidi == -321) return p[AliPID::kKaon]>0.5 && q==-1;
 143   else if (pidi ==  321) return p[AliPID::kKaon]>0.5 && q==1;
 144   else if (pidi == -2212) return p[AliPID::kProton]>0.5 && q==-1;
 145   else if (pidi ==  2212) return p[AliPID::kProton]>0.5 && q==1;
 146   else if (pidi == -11) return p[AliPID::kElectron]>0.5 && q==1;
 147   else if (pidi ==  11) return p[AliPID::kElectron]>0.5 && q==-1;
 148   else return 0;
 149 }
 150 //=============================================================================
 151 Bool_t AliUnicorEventAliceESD::PairGood(double p0, double the0, double phi0, double z0,
 152                                 double p1, double the1, double phi1, double z1) const {
 153
 154   // two-track separation cut
 155
 156   double dthe = the1-the0;
 157   if (fabs(dthe)>0.010) return kTRUE;
 158   double B = -0.5; // magnetic field in T, needed for helix; should be gotten in the proper way
 159   double pt0 = p0*sin(the0);
 160   double pt1 = p1*sin(the1);
 161   double r = 1.2; // we will calculate the distance between the two tracks at r=1.2 m
 162   // for (double r=0.8; r<2.5; r+=0.2) {
 163   double si0 = -0.3*B*z0*r/2/pt0;
 164   double si1 = -0.3*B*z1*r/2/pt1;
 165   if (fabs(si0)>=1.0) return kTRUE; // could be done better
 166   if (fabs(si1)>=1.0) return kTRUE;
 167   double dphi = phi1 - phi0 + asin(si1) - asin(si0);
 168   dphi = TVector2::Phi_mpi_pi(dphi);
 169   if (fabs(dphi)<0.020) return kFALSE;
 170   // }
 171   return kTRUE;
 172 }
 173 //=============================================================================