PWG2/SPECTRA/AliProtonAnalysis.cxx

   1 /**************************************************************************
   2  * Author: Panos Christakoglou.                                           *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 /* $Id$ */
  15
  16 //-----------------------------------------------------------------
  17 //                 AliProtonAnalysis class
  18 //   This is the class to deal with the proton analysis
  19 //   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
  20 //-----------------------------------------------------------------
  21 #include <Riostream.h>
  22 #include <TFile.h>
  23 #include <TSystem.h>
  24 #include <TH2F.h>
  25 #include <TH1D.h>
  26
  27 #include "AliProtonAnalysis.h"
  28
  29 #include <AliESDEvent.h>
  30 #include <AliLog.h>
  31
  32 ClassImp(AliProtonAnalysis)
  33
  34 //____________________________________________________________________//
  35 AliProtonAnalysis::AliProtonAnalysis() :
  36   TObject(),
  37   fNBinsY(0), fMinY(0), fMaxY(0),
  38   fNBinsPt(0), fMinPt(0), fMaxPt(0),
  39   fMinTPCClusters(0), fMinITSClusters(0),
  40   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  41   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  42   fMaxSigmaToVertex(0),
  43   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  44   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  45   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  46   fMaxSigmaToVertexFlag(kFALSE),
  47   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  48   fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  49   //Default constructor
  50   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
  51 }
  52
  53 //____________________________________________________________________//
  54 AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
  55   TObject(),
  56   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  57   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  58   fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  59   //Default constructor
  60
  61   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  62   fHistYPtProtons->SetStats(kTRUE);
  63   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  64   fHistYPtProtons->GetXaxis()->SetTitle("y");
  65   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
  66
  67   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  68   fHistYPtAntiProtons->SetStats(kTRUE);
  69   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  70   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
  71   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
  72 }
  73
  74 //____________________________________________________________________//
  75 AliProtonAnalysis::~AliProtonAnalysis() {
  76   //Default destructor
  77
  78 }
  79
  80 //____________________________________________________________________//
  81 void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
  82   fNBinsY = nbinsY;
  83   fMinY = fLowY;
  84   fMaxY = fHighY;
  85   fNBinsPt = nbinsPt;
  86   fMinPt = fLowPt;
  87   fMaxPt = fHighPt;
  88
  89   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  90   fHistYPtProtons->SetStats(kTRUE);
  91   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  92   fHistYPtProtons->GetXaxis()->SetTitle("y");
  93   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
  94
  95   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  96   fHistYPtAntiProtons->SetStats(kTRUE);
  97   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  98   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
  99   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 100 }
 101
 102 //____________________________________________________________________//
 103 void AliProtonAnalysis::ReadFromFile(const char* filename) {
 104   TFile *file = TFile::Open(filename);
 105   fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
 106   fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
 107   fHistYPtProtons->Sumw2();
 108   fHistYPtAntiProtons->Sumw2();
 109 }
 110
 111 //____________________________________________________________________//
 112 TH1D *AliProtonAnalysis::GetProtonYHistogram() {
 113   TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e");
 114   fYProtons->SetStats(kFALSE);
 115   fYProtons->GetYaxis()->SetTitle("dN/dy");
 116   fYProtons->SetTitle("dN/dy protons");
 117   fYProtons->SetMarkerStyle(kFullCircle);
 118   fYProtons->SetMarkerColor(4);
 119
 120   return fYProtons;
 121 }
 122
 123 //____________________________________________________________________//
 124 TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
 125   TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e");
 126   fYAntiProtons->SetStats(kFALSE);
 127   fYAntiProtons->GetYaxis()->SetTitle("dN/dy");
 128   fYAntiProtons->SetTitle("dN/dy antiprotons");
 129   fYAntiProtons->SetMarkerStyle(kFullCircle);
 130   fYAntiProtons->SetMarkerColor(4);
 131
 132   return fYAntiProtons;
 133 }
 134
 135 //____________________________________________________________________//
 136 TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
 137   TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
 138   fPtProtons->SetStats(kFALSE);
 139   fPtProtons->GetYaxis()->SetTitle("dN/dP_{T}");
 140   fPtProtons->SetTitle("dN/dPt protons");
 141   fPtProtons->SetMarkerStyle(kFullCircle);
 142   fPtProtons->SetMarkerColor(4);
 143
 144   return fPtProtons;
 145 }
 146
 147 //____________________________________________________________________//
 148 TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
 149   TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
 150   fPtAntiProtons->SetStats(kFALSE);
 151   fPtAntiProtons->GetYaxis()->SetTitle("dN/dP_{T}");
 152   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
 153   fPtAntiProtons->SetMarkerStyle(kFullCircle);
 154   fPtAntiProtons->SetMarkerColor(4);
 155
 156   return fPtAntiProtons;
 157 }
 158
 159 //____________________________________________________________________//
 160 TH1D *AliProtonAnalysis::GetYRatioHistogram() {
 161   TH1D *fYProtons = GetProtonYHistogram();
 162   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 163
 164   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 165   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
 166   hRatioY->SetMarkerStyle(kFullCircle);
 167   hRatioY->SetMarkerColor(4);
 168   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
 169   hRatioY->GetYaxis()->SetTitleOffset(1.4);
 170   hRatioY->GetXaxis()->SetTitle("y");
 171   hRatioY->GetXaxis()->SetTitleColor(1);
 172   hRatioY->SetStats(kFALSE);
 173
 174   return hRatioY;
 175 }
 176
 177 //____________________________________________________________________//
 178 TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
 179   TH1D *fPtProtons = GetProtonPtHistogram();
 180   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 181
 182   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 183   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
 184   hRatioPt->SetMarkerStyle(kFullCircle);
 185   hRatioPt->SetMarkerColor(4);
 186   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
 187   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
 188   hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 189   hRatioPt->GetXaxis()->SetTitleColor(1);
 190   hRatioPt->SetStats(kFALSE);
 191
 192   return hRatioPt;
 193 }
 194
 195 //____________________________________________________________________//
 196 TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
 197   TH1D *fYProtons = GetProtonYHistogram();
 198   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 199
 200   TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 201   hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
 202
 203   TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 204   hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
 205
 206   TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 207   hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
 208   hAsymmetryY->SetMarkerStyle(kFullCircle);
 209   hAsymmetryY->SetMarkerColor(4);
 210   hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
 211   hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
 212   hAsymmetryY->GetXaxis()->SetTitle("y");
 213   hAsymmetryY->GetXaxis()->SetTitleColor(1);
 214   hAsymmetryY->SetStats(kFALSE);
 215
 216   return hAsymmetryY;
 217 }
 218
 219 //____________________________________________________________________//
 220 TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
 221   TH1D *fPtProtons = GetProtonPtHistogram();
 222   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 223
 224   TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 225   hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
 226
 227   TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 228   hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
 229
 230   TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 231   hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
 232   hAsymmetryPt->SetMarkerStyle(kFullCircle);
 233   hAsymmetryPt->SetMarkerColor(4);
 234   hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
 235   hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
 236   hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 237   hAsymmetryPt->GetXaxis()->SetTitleColor(1);
 238   hAsymmetryPt->SetStats(kFALSE);
 239
 240   return hAsymmetryPt;
 241 }
 242
 243 //____________________________________________________________________//
 244 void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
 245   //Main analysis part
 246   Int_t nGoodTracks = fESD->GetNumberOfTracks();
 247   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 248     AliESDtrack* track = fESD->GetTrack(iTracks);
 249     if(IsAccepted(track)) {
 250       Double_t Pt = track->Pt();
 251
 252           //pid
 253       Double_t probability[5];
 254           track->GetESDpid(probability);
 255       Double_t rcc = 0.0;
 256       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += probability[i]*fPartFrac[i];
 257       if(rcc == 0.0) continue;
 258       Double_t w[5];
 259       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*fPartFrac[i]/rcc;
 260       Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 261       if(fParticleType == 4) {
 262         if(track->Charge() > 0) fHistYPtProtons->Fill(Rapidity(track),Pt);
 263         else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(Rapidity(track),Pt);
 264       }//proton check
 265     }//cuts
 266   }//track loop
 267 }
 268
 269 //____________________________________________________________________//
 270 Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
 271   // Checks if the track is excluded from the cuts
 272   Int_t  fIdxInt[200];
 273   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
 274   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
 275
 276   Float_t chi2PerClusterITS = -1;
 277   Float_t chi2PerClusterTPC = -1;
 278   if (nClustersTPC!=0)
 279     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 280
 281   Double_t extCov[15];
 282   track->GetExternalCovariance(extCov);
 283
 284   Double_t Pt = track->Pt();
 285   Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) +
 286                            TMath::Power(track->Py(),2) +
 287                            TMath::Power(track->Pz(),2));
 288
 289   if(fMinTPCClustersFlag)
 290     if(nClustersITS < fMinITSClusters) return kFALSE;
 291   if(fMinTPCClustersFlag)
 292     if(nClustersTPC < fMinTPCClusters) return kFALSE;
 293   if(fMaxChi2PerTPCClusterFlag)
 294     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
 295   if(fMaxChi2PerITSClusterFlag)
 296     if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
 297   if(fMaxCov11Flag)
 298     if(extCov[0] > fMaxCov11) return kFALSE;
 299   if(fMaxCov22Flag)
 300     if(extCov[2] > fMaxCov22) return kFALSE;
 301   if(fMaxCov33Flag)
 302     if(extCov[5] > fMaxCov33) return kFALSE;
 303   if(fMaxCov44Flag)
 304     if(extCov[9] > fMaxCov44) return kFALSE;
 305   if(fMaxCov55Flag)
 306     if(extCov[14] > fMaxCov55) return kFALSE;
 307   if(fMaxSigmaToVertexFlag)
 308     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
 309   if(fITSRefitFlag)
 310     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
 311   if(fTPCRefitFlag)
 312     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
 313
 314   if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
 315
 316   return kTRUE;
 317 }
 318
 319 //____________________________________________________________________//
 320 Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
 321   // Calculates the number of sigma to the vertex.
 322
 323   Float_t b[2];
 324   Float_t bRes[2];
 325   Float_t bCov[3];
 326   esdTrack->GetImpactParameters(b,bCov);
 327   if (bCov[0]<=0 || bCov[2]<=0) {
 328     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 329     bCov[0]=0; bCov[2]=0;
 330   }
 331   bRes[0] = TMath::Sqrt(bCov[0]);
 332   bRes[1] = TMath::Sqrt(bCov[2]);
 333
 334   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 335
 336   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 337
 338   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 339
 340   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 341
 342   return d;
 343 }
 344
 345 Double_t AliProtonAnalysis::Rapidity(AliESDtrack *track) {
 346   Double_t fMass = 9.38270000000000048e-01;
 347
 348   Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) +
 349                            TMath::Power(track->Py(),2) +
 350                                                    TMath::Power(track->Pz(),2));
 351   Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
 352   Double_t y = -999;
 353   if(energy != track->Pz())
 354     y = 0.5*TMath::Log((energy + track->Pz())/(energy - track->Pz()));
 355
 356   return y;
 357 }