PWG2/SPECTRA/AliProtonFeedDownAnalysis.cxx

   1 #include <Riostream.h>
   2 #include <TFile.h>
   3 #include <TSystem.h>
   4 #include <TF1.h>
   5 #include <TH2D.h>
   6 #include <TH1D.h>
   7 #include <TH1I.h>
   8 #include <TParticle.h>
   9
  10 #include <AliExternalTrackParam.h>
  11 #include <AliAODEvent.h>
  12 #include <AliESDEvent.h>
  13 //#include <AliLog.h>
  14 #include <AliPID.h>
  15 #include <AliStack.h>
  16 #include <AliCFContainer.h>
  17 #include <AliCFEffGrid.h>
  18 #include <AliCFDataGrid.h>
  19 //#include <AliESDVertex.h>
  20 class AliLog;
  21 class AliESDVertex;
  22
  23 #include "AliProtonFeedDownAnalysis.h"
  24 #include "AliProtonAnalysisBase.h"
  25
  26 ClassImp(AliProtonFeedDownAnalysis)
  27
  28 //____________________________________________________________________//
  29 AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis() :
  30   TObject(), fProtonAnalysisBase(0),
  31   fNBinsY(0), fMinY(0), fMaxY(0),
  32   fNBinsPt(0), fMinPt(0), fMaxPt(0),
  33   fProtonContainer(0), fAntiProtonContainer(0), fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
  34   {
  35   //Default constructor
  36  }
  37 //____________________________________________________________________//
  38 /*AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
  39   TObject(), fProtonAnalysisBase(0),
  40   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  41   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  42   fProtonContainer(0), fAntiProtonContainer(0),fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
  43   {
  44         //Default constructor
  45
  46         //setting up the containers
  47         Int_t iBin[2];
  48         iBin[0] = nbinsY;
  49         iBin[1] = nbinsPt;
  50         Double_t *binLimY = new Double_t[nbinsY+1];
  51         Double_t *binLimPt = new Double_t[nbinsPt+1];
  52         //values for bin lower bounds
  53         for(Int_t i = 0; i <= nbinsY; i++)
  54                 binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
  55         for(Int_t i = 0; i <= nbinsPt; i++)
  56                 binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
  57
  58         fProtonContainer = new AliCFContainer("containerProtons","container for protons",4,2,iBin);
  59         fProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
  60         fProtonContainer->SetBinLimits(1,binLimPt); //pT
  61         fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",4,2,iBin);
  62         fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
  63         fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
  64
  65
  66 }*/
  67 //____________________________________________________________________//
  68 AliProtonFeedDownAnalysis::~AliProtonFeedDownAnalysis()
  69 {
  70         //Default destructor
  71         if(fProtonAnalysisBase) delete fProtonAnalysisBase;
  72         if(fProtonContainer) delete fProtonContainer;
  73         if(fAntiProtonContainer) delete fAntiProtonContainer;
  74         if(fweightfunction) delete fweightfunction;
  75         if(fLambda) delete fLambda;
  76         if(fLambdaweighted) delete fLambdaweighted;
  77         if(fAntiLambda) delete fAntiLambda;
  78         if(fAntiLambdaweighted) delete fAntiLambdaweighted;
  79 }
  80 //____________________________________________________________________//
  81 void AliProtonFeedDownAnalysis::InitAnalysisHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt)
  82 {
  83         //Initializes the histograms
  84         fNBinsY = nbinsY;
  85         fMinY = fLowY;
  86         fMaxY = fHighY;
  87         fNBinsPt = nbinsPt;
  88         fMinPt = fLowPt;
  89         fMaxPt = fHighPt;
  90
  91
  92         //setting up the containers
  93         Int_t iBin[2];
  94         iBin[0] = nbinsY;
  95         iBin[1] = nbinsPt;
  96         Double_t *binLimY = new Double_t[nbinsY+1];
  97         Double_t *binLimPt = new Double_t[nbinsPt+1];
  98         //values for bin lower bounds
  99         for(Int_t i = 0; i <= nbinsY; i++)
 100                 binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
 101         for(Int_t i = 0; i <= nbinsPt; i++)
 102                 binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
 103
 104         fProtonContainer = new AliCFContainer("containerProtons","container for protons",kNSteps,2,iBin);
 105         fProtonContainer->SetBinLimits(0,binLimY); //rapidity
 106         fProtonContainer->SetBinLimits(1,binLimPt); //pT
 107         fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",kNSteps,2,iBin);
 108         fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
 109         fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
 110         /*fLambda=new TH1F("Lambda","Lambda",35,0.5,4.0);
 111         fLambdaweighted=new TH1F("Lambdaweighted","Lambdaweighted",35,0.5,4.0);
 112         fAntiLambda=new TH1F("AntiLambda","AntiLambda",35,0.5,4.0);
 113         fAntiLambdaweighted=new TH1F("AntiLambdaweighted","AntiLambdaweighted",35,0.5,4.0);*/
 114         fLambda=new TH2F("Lambda","Lambda",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 115         fLambdaweighted=new TH2F("Lambdaweighted","Lambdaweighted",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 116         fAntiLambda=new TH2F("AntiLambda","AntiLambda",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 117         fAntiLambdaweighted=new TH2F("AntiLambdaweighted","AntiLambdaweighted",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 118         fLambda->GetYaxis()->SetTitle("P_{T} [GeV/c]");
 119         fLambdaweighted->GetYaxis()->SetTitle("P_{T} [GeV/c]");
 120         fAntiLambda->GetYaxis()->SetTitle("P_{T} [GeV/c]");
 121         fAntiLambdaweighted->GetYaxis()->SetTitle("P_{T} [GeV/c]");
 122         if(fProtonAnalysisBase->GetEtaMode())
 123         {
 124                 fLambda->GetXaxis()->SetTitle("#eta");
 125                 fLambdaweighted->GetXaxis()->SetTitle("#eta");
 126                 fAntiLambda->GetXaxis()->SetTitle("#eta");
 127                 fAntiLambdaweighted->GetXaxis()->SetTitle("#eta");
 128         }
 129         else
 130         {
 131                 fLambda->GetXaxis()->SetTitle("y");
 132                 fLambdaweighted->GetXaxis()->SetTitle("y");
 133                 fAntiLambda->GetXaxis()->SetTitle("y");
 134                 fAntiLambdaweighted->GetXaxis()->SetTitle("y");
 135         }
 136 }
 137 //_________________________________________________________________________//
 138 void AliProtonFeedDownAnalysis::Analyze(AliESDEvent *esd, const AliESDVertex *vertex,AliStack *stack)
 139 {
 140         Int_t nTracks = 0;
 141         Int_t nIdentifiedProtons = 0, nIdentifiedAntiProtons = 0;
 142         Int_t nSurvivedProtons = 0, nSurvivedAntiProtons = 0;
 143
 144         //fHistEvents->Fill(0); //number of analyzed events
 145         Double_t containerInput[2] ;
 146         Double_t gPt = 0.0, gP = 0.0;
 147         nTracks = esd->GetNumberOfTracks();
 148         Float_t weight;
 149         Int_t trackflag;
 150         for(Int_t iTracks = 0; iTracks < nTracks; iTracks++)
 151         {
 152
 153                 AliESDtrack* track = esd->GetTrack(iTracks);
 154                 AliESDtrack trackTPC;
 155                 Int_t label= track->GetLabel();
 156         /*      Int_t trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
 157                 if (trackflag!=0)
 158                         weight=GetWeightforProton(label,stack);
 159                 else
 160                         weight=1.0;     */
 161
 162                 if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid))
 163                 {
 164                         AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 165                         if(!tpcTrack)
 166                                  continue;
 167                         gPt = tpcTrack->Pt();
 168                         gP = tpcTrack->P();
 169                         if(fProtonAnalysisBase->IsProton(track))
 170                         {
 171                                 trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
 172                                 if (trackflag!=0)
 173                                         weight=GetWeightforProton(label,stack);
 174                                 else
 175                                         weight=1.0;
 176                                 if(tpcTrack->Charge() > 0)
 177                                 {
 178                                         nIdentifiedProtons += 1;
 179                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track))
 180                                                 continue;//track cuts
 181                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 182                                                 continue; //track outside the analyzed y-Pt
 183                                         nSurvivedProtons += 1;
 184                                         if(fProtonAnalysisBase->GetEtaMode())
 185                                         {
 186                                                 containerInput[0] = tpcTrack->Eta();
 187                                         }
 188                                         else
 189                                         {
 190                                                 //fill the container
 191                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
 192                                         }
 193                                         containerInput[1] = gPt;
 194                                         fProtonContainer->Fill(containerInput,kSelected,weight);
 195                                         //Feed-down check
 196                                         if (trackflag==1)
 197                                         {
 198                                                 fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 199                                                 TParticle *particle  = stack->Particle(label);
 200                                                 Int_t lmother=particle->GetFirstMother();
 201                                                 TParticle *mparticle=stack->Particle(lmother);
 202                                                 Double_t ptlambda= mparticle->Pt();
 203                                                 Double_t eta_or_y=0.0;
 204                                                 if(fProtonAnalysisBase->GetEtaMode())
 205                                                         eta_or_y=mparticle->Eta();
 206                                                 else
 207                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 208                                                 fLambda->Fill(eta_or_y, ptlambda);
 209                                                 fLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 210                                         }
 211
 212                                 }//protons
 213                                 else if(tpcTrack->Charge() < 0)
 214                                 {
 215                                         nIdentifiedAntiProtons += 1;
 216                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track))
 217                                                 continue;//track cuts
 218                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 219                                                 continue; //track outside the analyzed y-Pt
 220                                         nSurvivedAntiProtons += 1;
 221                                         if(fProtonAnalysisBase->GetEtaMode())
 222                                         {
 223                                                 containerInput[0] = tpcTrack->Eta();
 224                                         }
 225                                         else
 226                                         {
 227                                                 //fill the container
 228                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
 229                                         }
 230                                         containerInput[1] = gPt;
 231                                         fAntiProtonContainer->Fill(containerInput,kSelected,weight);
 232                                         //Feed-down check
 233                                         if(trackflag==-1)
 234                                         {
 235                                                 fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 236                                                 TParticle *particle  = stack->Particle(label);
 237                                                 Int_t lmother=particle->GetFirstMother();
 238                                                 TParticle *mparticle=stack->Particle(lmother);
 239                                                 Double_t ptlambda= mparticle->Pt();
 240                                                 Double_t eta_or_y=0.0;
 241                                                 if(fProtonAnalysisBase->GetEtaMode())
 242                                                         eta_or_y=mparticle->Eta();
 243                                                 else
 244                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 245                                                 fAntiLambda->Fill(eta_or_y, ptlambda);
 246                                                 fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 247                                         }
 248                                 }//antiprotons
 249                         }//proton check
 250                 }//TPC only tracks
 251                 else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal)
 252                 {
 253                         gPt = track->Pt();
 254                         gP = track->P();
 255                         if(fProtonAnalysisBase->IsProton(track))
 256                         {
 257                                 trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
 258                                 if (trackflag!=0)
 259                                         weight=GetWeightforProton(label,stack);
 260                                 else
 261                                         weight=1.0;
 262                                 if(track->Charge() > 0)
 263                                 {
 264                                         nIdentifiedProtons += 1;
 265                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track))
 266                                                 continue;//track cuts
 267                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 268                                                 continue; //track outside the analyzed y-Pt
 269                                         nSurvivedProtons += 1;
 270                                         if(fProtonAnalysisBase->GetEtaMode())
 271                                         {
 272                                                 containerInput[0] = track->Eta();
 273                                         }
 274                                         else
 275                                         {
 276                                                 //fill the container
 277                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
 278                                         }
 279                                         containerInput[1] = gPt;
 280                                         fProtonContainer->Fill(containerInput,kSelected,weight);
 281                                          //Feed-down check
 282                                         if (trackflag==1)
 283                                         {
 284                                                 fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 285                                                 TParticle *particle  = stack->Particle(label);
 286                                                 Int_t lmother=particle->GetFirstMother();
 287                                                 TParticle *mparticle=stack->Particle(lmother);
 288                                                 Double_t ptlambda= mparticle->Pt();
 289                                                 Double_t eta_or_y=0.0;
 290                                                 if(fProtonAnalysisBase->GetEtaMode())
 291                                                         eta_or_y=mparticle->Eta();
 292                                                 else
 293                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 294                                                 fLambda->Fill(eta_or_y, ptlambda);
 295                                                 fLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 296                                         }
 297                                 }//protons
 298                                 else if(track->Charge() < 0)
 299                                 {
 300                                         nIdentifiedAntiProtons += 1;
 301                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track))
 302                                                 continue;//track cuts
 303                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 304                                                 continue; //track outside the analyzed y-Pt
 305                                         nSurvivedAntiProtons += 1;
 306                                         if(fProtonAnalysisBase->GetEtaMode())
 307                                         {
 308                                                 containerInput[0] = track->Eta();
 309                                         }
 310                                         else
 311                                         {
 312                                                 //fill the container
 313                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
 314                                         }
 315                                         containerInput[1] = gPt;
 316                                         fAntiProtonContainer->Fill(containerInput,kSelected,weight);
 317                                         if(trackflag==-1)
 318                                         {
 319                                                 fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 320                                                 TParticle *particle  = stack->Particle(label);
 321                                                 Int_t lmother=particle->GetFirstMother();
 322                                                 TParticle *mparticle=stack->Particle(lmother);
 323                                                 Double_t ptlambda= mparticle->Pt();
 324                                                 Double_t eta_or_y=0.0;
 325                                                 if(fProtonAnalysisBase->GetEtaMode())
 326                                                         eta_or_y=mparticle->Eta();
 327                                                 else
 328                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 329                                                 fAntiLambda->Fill(eta_or_y, ptlambda);
 330                                                 fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 331                                         }
 332                                 }//antiprotons
 333                         }//proton check
 334                 }//combined tracking
 335         }//track loop
 336
 337         if(fProtonAnalysisBase->GetDebugMode())
 338                 Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
 339
 340 }
 341 //_________________________________________________________________________//
 342 void AliProtonFeedDownAnalysis::Analyze(AliAODEvent *fAOD)
 343 {
 344   // Analysis from AOD: to be implemented...
 345   // in the near future.
 346   fAOD->Print();
 347
 348 }
 349 //_________________________________________________________________________//
 350 void AliProtonFeedDownAnalysis::Analyze(AliStack *stack)
 351 {
 352          Double_t containerInput[2] ;
 353          Float_t weight;
 354         for (Int_t ipart=0; ipart<stack->GetNtrack(); ipart++)
 355         {
 356                 TParticle *particle  = stack->Particle(ipart);
 357                 Int_t code=particle->GetPdgCode();
 358                 /*if (code==3122)
 359                 {
 360                         fLambda->Fill(particle->Pt());
 361                         fLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));
 362                 }
 363                 if (code==-3122)
 364                 {
 365                         fAntiLambda->Fill(particle->Pt());
 366                         fAntiLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));
 367                 }*/
 368                 if (TMath::Abs(code)==2212)
 369                 {
 370                         Int_t trackflag=LambdaIsMother(ipart,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
 371                         if (trackflag!=0)
 372                                 weight=GetWeightforProton(ipart,stack);
 373                         else
 374                                 weight=1.0;
 375                         if(particle->GetUniqueID() == 13) //recject hadronic inter.
 376                                 continue;
 377                         if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt))
 378                                 continue;
 379                         if(fProtonAnalysisBase->GetEtaMode())
 380                         {
 381                                 if((particle->Eta()> fMaxY)||(particle->Eta()< fMinY))
 382                                         continue;
 383                         }
 384                         else
 385                         {
 386                                 if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz())< fMinY))
 387                                         continue;
 388                         }
 389                         if(fProtonAnalysisBase->GetEtaMode())
 390                         {
 391                                 containerInput[0] =particle->Eta();
 392                         }
 393                         else
 394                         {
 395                                 containerInput[0] = fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz());
 396                         }
 397                         //containerInput[0] = particle->Eta() ;
 398                         containerInput[1] = (Float_t)particle->Pt();
 399                         if (particle->GetPdgCode()==2212)
 400                         {
 401                                 fProtonContainer->Fill(containerInput,kAll,weight);
 402                                 if(ipart<stack->GetNprimary())
 403                                 {
 404                                         fProtonContainer->Fill(containerInput,kPrimary,weight);
 405                                 }
 406                                 else
 407                                 {
 408                                         if(trackflag==1)
 409                                                 fProtonContainer->Fill(containerInput,kFromLambda,weight);
 410                                 }
 411                         }
 412                         if (particle->GetPdgCode()==-2212)
 413                         {
 414                                 fAntiProtonContainer->Fill(containerInput,kAll,weight);
 415                                 if(ipart<stack->GetNprimary())
 416                                 {
 417                                         fAntiProtonContainer->Fill(containerInput,kPrimary,weight);
 418                                 }
 419                                 else
 420                                 {
 421                                         if(trackflag==-1)
 422                                                 fAntiProtonContainer->Fill(containerInput,kFromLambda,weight);
 423                                 }
 424                         }
 425                 }
 426         }
 427 }
 428 //______________________________________________________________________________________________
 429 Int_t AliProtonFeedDownAnalysis::LambdaIsMother(Int_t number, AliStack *stack)
 430 {
 431          if(number<0)
 432                 return 0;
 433         TParticle *particle  = stack->Particle(number);
 434         Int_t motherPDG=0;
 435         Int_t lmother=-1;
 436         lmother=particle->GetFirstMother();
 437         if (lmother<0)
 438                 return 0;
 439         TParticle *mparticle=stack->Particle(lmother);
 440         motherPDG=mparticle->GetPdgCode();
 441         if(motherPDG==3122)
 442                 return 1;
 443         if(motherPDG==-3122)
 444                 return -1;
 445         return 0;
 446 }
 447 //___________________________________________________________________________________________
 448 Float_t AliProtonFeedDownAnalysis::GetWeightforProton(Int_t number,AliStack *stack)
 449 {
 450          if(number<0)
 451                 return 1.0;
 452         TParticle *particle  = stack->Particle(number);
 453         Int_t lmother=particle->GetFirstMother();
 454         TParticle *mparticle=stack->Particle(lmother);
 455         return  fweightfunction->Eval(mparticle->Pt());
 456 }
 457 //______________________________________________________________________________________________
 458 Float_t AliProtonFeedDownAnalysis::GetWeightforLambda(Float_t pt)
 459 {
 460         return  fweightfunction->Eval(pt);
 461 }