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 <TF1.h>
  25 #include <TH2F.h>
  26 #include <TH1D.h>
  27
  28 #include "AliProtonAnalysis.h"
  29
  30 #include <AliAODEvent.h>
  31 #include <AliESDEvent.h>
  32 #include <AliLog.h>
  33 #include <AliPID.h>
  34
  35 ClassImp(AliProtonAnalysis)
  36
  37 //____________________________________________________________________//
  38 AliProtonAnalysis::AliProtonAnalysis() :
  39   TObject(),
  40   fNBinsY(0), fMinY(0), fMaxY(0),
  41   fNBinsPt(0), fMinPt(0), fMaxPt(0),
  42   fMinTPCClusters(0), fMinITSClusters(0),
  43   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  44   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  45   fMaxSigmaToVertex(0),
  46   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  47   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  48   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  49   fMaxSigmaToVertexFlag(kFALSE),
  50   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  51   fFunctionProbabilityFlag(kFALSE),
  52   fElectronFunction(0), fMuonFunction(0),
  53   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  54   fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  55   //Default constructor
  56   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
  57 }
  58
  59 //____________________________________________________________________//
  60 AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
  61   TObject(),
  62   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  63   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  64   fMinTPCClusters(0), fMinITSClusters(0),
  65   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  66   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  67   fMaxSigmaToVertex(0),
  68   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  69   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  70   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  71   fMaxSigmaToVertexFlag(kFALSE),
  72   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  73   fFunctionProbabilityFlag(kFALSE),
  74   fElectronFunction(0), fMuonFunction(0),
  75   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  76   fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  77   //Default constructor
  78
  79   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  80   fHistYPtProtons->SetStats(kTRUE);
  81   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  82   fHistYPtProtons->GetXaxis()->SetTitle("y");
  83   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
  84
  85   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  86   fHistYPtAntiProtons->SetStats(kTRUE);
  87   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  88   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
  89   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
  90 }
  91
  92 //____________________________________________________________________//
  93 AliProtonAnalysis::~AliProtonAnalysis() {
  94   //Default destructor
  95
  96 }
  97
  98 //____________________________________________________________________//
  99 void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
 100   fNBinsY = nbinsY;
 101   fMinY = fLowY;
 102   fMaxY = fHighY;
 103   fNBinsPt = nbinsPt;
 104   fMinPt = fLowPt;
 105   fMaxPt = fHighPt;
 106
 107   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 108   fHistYPtProtons->SetStats(kTRUE);
 109   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 110   fHistYPtProtons->GetXaxis()->SetTitle("y");
 111   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
 112
 113   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 114   fHistYPtAntiProtons->SetStats(kTRUE);
 115   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 116   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
 117   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 118 }
 119
 120 //____________________________________________________________________//
 121 void AliProtonAnalysis::ReadFromFile(const char* filename) {
 122   TFile *file = TFile::Open(filename);
 123   fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
 124   fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
 125   fHistYPtProtons->Sumw2();
 126   fHistYPtAntiProtons->Sumw2();
 127 }
 128
 129 //____________________________________________________________________//
 130 TH1D *AliProtonAnalysis::GetProtonYHistogram() {
 131   TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e");
 132   fYProtons->SetStats(kFALSE);
 133   fYProtons->GetYaxis()->SetTitle("dN/dy");
 134   fYProtons->SetTitle("dN/dy protons");
 135   fYProtons->SetMarkerStyle(kFullCircle);
 136   fYProtons->SetMarkerColor(4);
 137
 138   return fYProtons;
 139 }
 140
 141 //____________________________________________________________________//
 142 TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
 143   TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e");
 144   fYAntiProtons->SetStats(kFALSE);
 145   fYAntiProtons->GetYaxis()->SetTitle("dN/dy");
 146   fYAntiProtons->SetTitle("dN/dy antiprotons");
 147   fYAntiProtons->SetMarkerStyle(kFullCircle);
 148   fYAntiProtons->SetMarkerColor(4);
 149
 150   return fYAntiProtons;
 151 }
 152
 153 //____________________________________________________________________//
 154 TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
 155   TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
 156   fPtProtons->SetStats(kFALSE);
 157   fPtProtons->GetYaxis()->SetTitle("dN/dP_{T}");
 158   fPtProtons->SetTitle("dN/dPt protons");
 159   fPtProtons->SetMarkerStyle(kFullCircle);
 160   fPtProtons->SetMarkerColor(4);
 161
 162   return fPtProtons;
 163 }
 164
 165 //____________________________________________________________________//
 166 TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
 167   TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
 168   fPtAntiProtons->SetStats(kFALSE);
 169   fPtAntiProtons->GetYaxis()->SetTitle("dN/dP_{T}");
 170   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
 171   fPtAntiProtons->SetMarkerStyle(kFullCircle);
 172   fPtAntiProtons->SetMarkerColor(4);
 173
 174   return fPtAntiProtons;
 175 }
 176
 177 //____________________________________________________________________//
 178 TH1D *AliProtonAnalysis::GetYRatioHistogram() {
 179   TH1D *fYProtons = GetProtonYHistogram();
 180   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 181
 182   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 183   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
 184   hRatioY->SetMarkerStyle(kFullCircle);
 185   hRatioY->SetMarkerColor(4);
 186   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
 187   hRatioY->GetYaxis()->SetTitleOffset(1.4);
 188   hRatioY->GetXaxis()->SetTitle("y");
 189   hRatioY->GetXaxis()->SetTitleColor(1);
 190   hRatioY->SetStats(kFALSE);
 191
 192   return hRatioY;
 193 }
 194
 195 //____________________________________________________________________//
 196 TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
 197   TH1D *fPtProtons = GetProtonPtHistogram();
 198   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 199
 200   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 201   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
 202   hRatioPt->SetMarkerStyle(kFullCircle);
 203   hRatioPt->SetMarkerColor(4);
 204   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
 205   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
 206   hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 207   hRatioPt->GetXaxis()->SetTitleColor(1);
 208   hRatioPt->SetStats(kFALSE);
 209
 210   return hRatioPt;
 211 }
 212
 213 //____________________________________________________________________//
 214 TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
 215   TH1D *fYProtons = GetProtonYHistogram();
 216   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 217
 218   TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 219   hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
 220
 221   TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 222   hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
 223
 224   TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 225   hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
 226   hAsymmetryY->SetMarkerStyle(kFullCircle);
 227   hAsymmetryY->SetMarkerColor(4);
 228   hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
 229   hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
 230   hAsymmetryY->GetXaxis()->SetTitle("y");
 231   hAsymmetryY->GetXaxis()->SetTitleColor(1);
 232   hAsymmetryY->SetStats(kFALSE);
 233
 234   return hAsymmetryY;
 235 }
 236
 237 //____________________________________________________________________//
 238 TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
 239   TH1D *fPtProtons = GetProtonPtHistogram();
 240   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 241
 242   TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 243   hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
 244
 245   TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 246   hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
 247
 248   TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 249   hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
 250   hAsymmetryPt->SetMarkerStyle(kFullCircle);
 251   hAsymmetryPt->SetMarkerColor(4);
 252   hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
 253   hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
 254   hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 255   hAsymmetryPt->GetXaxis()->SetTitleColor(1);
 256   hAsymmetryPt->SetStats(kFALSE);
 257
 258   return hAsymmetryPt;
 259 }
 260
 261 //____________________________________________________________________//
 262 Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
 263   Double_t partFrac=0;
 264   if(fFunctionProbabilityFlag) {
 265     if(i == 0) partFrac = fElectronFunction->Eval(p);
 266     if(i == 1) partFrac = fMuonFunction->Eval(p);
 267     if(i == 2) partFrac = fPionFunction->Eval(p);
 268     if(i == 3) partFrac = fKaonFunction->Eval(p);
 269     if(i == 4) partFrac = fProtonFunction->Eval(p);
 270   }
 271   else partFrac = fPartFrac[i];
 272
 273   return partFrac;
 274 }
 275
 276 //____________________________________________________________________//
 277 void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
 278   //Main analysis part
 279   Int_t nGoodTracks = fESD->GetNumberOfTracks();
 280   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 281     AliESDtrack* track = fESD->GetTrack(iTracks);
 282     if(IsAccepted(track)) {
 283       Double_t Pt = track->Pt();
 284       Double_t P = track->P();
 285
 286       //pid
 287       Double_t probability[5];
 288           track->GetESDpid(probability);
 289       Double_t rcc = 0.0;
 290       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += probability[i]*GetParticleFraction(i,P);
 291       if(rcc == 0.0) continue;
 292       Double_t w[5];
 293       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
 294       Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 295       if(fParticleType == 4) {
 296         if(track->Charge() > 0) fHistYPtProtons->Fill(Rapidity(track),Pt);
 297         else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(Rapidity(track),Pt);
 298       }//proton check
 299     }//cuts
 300   }//track loop
 301 }
 302
 303 //____________________________________________________________________//
 304 void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
 305   //Main analysis part
 306   Int_t nGoodTracks = fAOD->GetNumberOfTracks();
 307   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 308     AliAODTrack* track = fAOD->GetTrack(iTracks);
 309     Double_t Pt = track->Pt();
 310     Double_t P = track->P();
 311
 312     //pid
 313     Double_t probability[10];
 314     track->GetPID(probability);
 315     Double_t rcc = 0.0;
 316     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
 317     if(rcc == 0.0) continue;
 318     Double_t w[10];
 319     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
 320     Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
 321     if(fParticleType == 4) {
 322       if(track->Charge() > 0) fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
 323       else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
 324     }//proton check
 325   }//track loop
 326 }
 327
 328 //____________________________________________________________________//
 329 Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
 330   // Checks if the track is excluded from the cuts
 331   Int_t  fIdxInt[200];
 332   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
 333   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
 334
 335   Float_t chi2PerClusterITS = -1;
 336   Float_t chi2PerClusterTPC = -1;
 337   if (nClustersTPC!=0)
 338     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 339
 340   Double_t extCov[15];
 341   track->GetExternalCovariance(extCov);
 342
 343   Double_t Pt = track->Pt();
 344
 345   if(fMinITSClustersFlag)
 346     if(nClustersITS < fMinITSClusters) return kFALSE;
 347   if(fMinTPCClustersFlag)
 348     if(nClustersTPC < fMinTPCClusters) return kFALSE;
 349   if(fMaxChi2PerTPCClusterFlag)
 350     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
 351   if(fMaxChi2PerITSClusterFlag)
 352     if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
 353   if(fMaxCov11Flag)
 354     if(extCov[0] > fMaxCov11) return kFALSE;
 355   if(fMaxCov22Flag)
 356     if(extCov[2] > fMaxCov22) return kFALSE;
 357   if(fMaxCov33Flag)
 358     if(extCov[5] > fMaxCov33) return kFALSE;
 359   if(fMaxCov44Flag)
 360     if(extCov[9] > fMaxCov44) return kFALSE;
 361   if(fMaxCov55Flag)
 362     if(extCov[14] > fMaxCov55) return kFALSE;
 363   if(fMaxSigmaToVertexFlag)
 364     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
 365   if(fITSRefitFlag)
 366     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
 367   if(fTPCRefitFlag)
 368     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
 369
 370   if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
 371
 372   return kTRUE;
 373 }
 374
 375 //____________________________________________________________________//
 376 Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
 377   // Calculates the number of sigma to the vertex.
 378
 379   Float_t b[2];
 380   Float_t bRes[2];
 381   Float_t bCov[3];
 382   esdTrack->GetImpactParameters(b,bCov);
 383   if (bCov[0]<=0 || bCov[2]<=0) {
 384     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 385     bCov[0]=0; bCov[2]=0;
 386   }
 387   bRes[0] = TMath::Sqrt(bCov[0]);
 388   bRes[1] = TMath::Sqrt(bCov[2]);
 389
 390   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 391
 392   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 393
 394   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 395
 396   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 397
 398   return d;
 399 }
 400
 401 Double_t AliProtonAnalysis::Rapidity(AliESDtrack *track) {
 402   Double_t fMass = 9.38270000000000048e-01;
 403
 404   Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) +
 405                            TMath::Power(track->Py(),2) +
 406                                                    TMath::Power(track->Pz(),2));
 407   Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
 408   Double_t y = -999;
 409   if(energy != track->Pz())
 410     y = 0.5*TMath::Log((energy + track->Pz())/(energy - track->Pz()));
 411
 412   return y;
 413 }