PWGLF/SPECTRA/AntiprotonToProton/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->IsPrimary(esd,vertex,track))
 180                                           continue;
 181                                         if(!fProtonAnalysisBase->IsAccepted(track))
 182                                                 continue;//track cuts
 183                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 184                                                 continue; //track outside the analyzed y-Pt
 185                                         nSurvivedProtons += 1;
 186                                         if(fProtonAnalysisBase->GetEtaMode())
 187                                         {
 188                                                 containerInput[0] = tpcTrack->Eta();
 189                                         }
 190                                         else
 191                                         {
 192                                                 //fill the container
 193                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
 194                                         }
 195                                         containerInput[1] = gPt;
 196                                         fProtonContainer->Fill(containerInput,kSelected,weight);
 197                                         //Feed-down check
 198                                         if (trackflag==1)
 199                                         {
 200                                                 fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 201                                                 TParticle *particle  = stack->Particle(label);
 202                                                 Int_t lmother=particle->GetFirstMother();
 203                                                 TParticle *mparticle=stack->Particle(lmother);
 204                                                 Double_t ptlambda= mparticle->Pt();
 205                                                 Double_t eta_or_y=0.0;
 206                                                 if(fProtonAnalysisBase->GetEtaMode())
 207                                                         eta_or_y=mparticle->Eta();
 208                                                 else
 209                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 210                                                 fLambda->Fill(eta_or_y, ptlambda);
 211                                                 fLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 212                                         }
 213
 214                                 }//protons
 215                                 else if(tpcTrack->Charge() < 0)
 216                                 {
 217                                         nIdentifiedAntiProtons += 1;
 218                                         if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track))
 219                                           continue;
 220                                         if(!fProtonAnalysisBase->IsAccepted(track))
 221                                                 continue;//track cuts
 222                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 223                                                 continue; //track outside the analyzed y-Pt
 224                                         nSurvivedAntiProtons += 1;
 225                                         if(fProtonAnalysisBase->GetEtaMode())
 226                                         {
 227                                                 containerInput[0] = tpcTrack->Eta();
 228                                         }
 229                                         else
 230                                         {
 231                                                 //fill the container
 232                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
 233                                         }
 234                                         containerInput[1] = gPt;
 235                                         fAntiProtonContainer->Fill(containerInput,kSelected,weight);
 236                                         //Feed-down check
 237                                         if(trackflag==-1)
 238                                         {
 239                                                 fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 240                                                 TParticle *particle  = stack->Particle(label);
 241                                                 Int_t lmother=particle->GetFirstMother();
 242                                                 TParticle *mparticle=stack->Particle(lmother);
 243                                                 Double_t ptlambda= mparticle->Pt();
 244                                                 Double_t eta_or_y=0.0;
 245                                                 if(fProtonAnalysisBase->GetEtaMode())
 246                                                         eta_or_y=mparticle->Eta();
 247                                                 else
 248                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 249                                                 fAntiLambda->Fill(eta_or_y, ptlambda);
 250                                                 fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 251                                         }
 252                                 }//antiprotons
 253                         }//proton check
 254                 }//TPC only tracks
 255                 else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal)
 256                 {
 257                         gPt = track->Pt();
 258                         gP = track->P();
 259                         if(fProtonAnalysisBase->IsProton(track))
 260                         {
 261                                 trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
 262                                 if (trackflag!=0)
 263                                         weight=GetWeightforProton(label,stack);
 264                                 else
 265                                         weight=1.0;
 266                                 if(track->Charge() > 0)
 267                                 {
 268                                         nIdentifiedProtons += 1;
 269                                         if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track))
 270                                           continue;
 271                                         if(!fProtonAnalysisBase->IsAccepted(track))
 272                                                 continue;//track cuts
 273                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 274                                                 continue; //track outside the analyzed y-Pt
 275                                         nSurvivedProtons += 1;
 276                                         if(fProtonAnalysisBase->GetEtaMode())
 277                                         {
 278                                                 containerInput[0] = track->Eta();
 279                                         }
 280                                         else
 281                                         {
 282                                                 //fill the container
 283                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
 284                                         }
 285                                         containerInput[1] = gPt;
 286                                         fProtonContainer->Fill(containerInput,kSelected,weight);
 287                                          //Feed-down check
 288                                         if (trackflag==1)
 289                                         {
 290                                                 fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 291                                                 TParticle *particle  = stack->Particle(label);
 292                                                 Int_t lmother=particle->GetFirstMother();
 293                                                 TParticle *mparticle=stack->Particle(lmother);
 294                                                 Double_t ptlambda= mparticle->Pt();
 295                                                 Double_t eta_or_y=0.0;
 296                                                 if(fProtonAnalysisBase->GetEtaMode())
 297                                                         eta_or_y=mparticle->Eta();
 298                                                 else
 299                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 300                                                 fLambda->Fill(eta_or_y, ptlambda);
 301                                                 fLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 302                                         }
 303                                 }//protons
 304                                 else if(track->Charge() < 0)
 305                                 {
 306                                         nIdentifiedAntiProtons += 1;
 307                                         if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track))
 308                                           continue;
 309                                         if(!fProtonAnalysisBase->IsAccepted(track))
 310                                                 continue;//track cuts
 311                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
 312                                                 continue; //track outside the analyzed y-Pt
 313                                         nSurvivedAntiProtons += 1;
 314                                         if(fProtonAnalysisBase->GetEtaMode())
 315                                         {
 316                                                 containerInput[0] = track->Eta();
 317                                         }
 318                                         else
 319                                         {
 320                                                 //fill the container
 321                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
 322                                         }
 323                                         containerInput[1] = gPt;
 324                                         fAntiProtonContainer->Fill(containerInput,kSelected,weight);
 325                                         if(trackflag==-1)
 326                                         {
 327                                                 fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
 328                                                 TParticle *particle  = stack->Particle(label);
 329                                                 Int_t lmother=particle->GetFirstMother();
 330                                                 TParticle *mparticle=stack->Particle(lmother);
 331                                                 Double_t ptlambda= mparticle->Pt();
 332                                                 Double_t eta_or_y=0.0;
 333                                                 if(fProtonAnalysisBase->GetEtaMode())
 334                                                         eta_or_y=mparticle->Eta();
 335                                                 else
 336                                                         eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
 337                                                 fAntiLambda->Fill(eta_or_y, ptlambda);
 338                                                 fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);
 339                                         }
 340                                 }//antiprotons
 341                         }//proton check
 342                 }//combined tracking
 343         }//track loop
 344
 345         if(fProtonAnalysisBase->GetDebugMode())
 346                 Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
 347
 348 }
 349 //_________________________________________________________________________//
 350 void AliProtonFeedDownAnalysis::Analyze(AliAODEvent *fAOD)
 351 {
 352   // Analysis from AOD: to be implemented...
 353   // in the near future.
 354   fAOD->Print();
 355
 356 }
 357 //_________________________________________________________________________//
 358 void AliProtonFeedDownAnalysis::Analyze(AliStack *stack)
 359 {
 360          Double_t containerInput[2] ;
 361          Float_t weight;
 362         for (Int_t ipart=0; ipart<stack->GetNtrack(); ipart++)
 363         {
 364                 TParticle *particle  = stack->Particle(ipart);
 365                 Int_t code=particle->GetPdgCode();
 366                 /*if (code==3122)
 367                 {
 368                         fLambda->Fill(particle->Pt());
 369                         fLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));
 370                 }
 371                 if (code==-3122)
 372                 {
 373                         fAntiLambda->Fill(particle->Pt());
 374                         fAntiLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));
 375                 }*/
 376                 if (TMath::Abs(code)==2212)
 377                 {
 378                         Int_t trackflag=LambdaIsMother(ipart,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
 379                         if (trackflag!=0)
 380                                 weight=GetWeightforProton(ipart,stack);
 381                         else
 382                                 weight=1.0;
 383                         if(particle->GetUniqueID() == 13) //recject hadronic inter.
 384                                 continue;
 385                         if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt))
 386                                 continue;
 387                         if(fProtonAnalysisBase->GetEtaMode())
 388                         {
 389                                 if((particle->Eta()> fMaxY)||(particle->Eta()< fMinY))
 390                                         continue;
 391                         }
 392                         else
 393                         {
 394                                 if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz())< fMinY))
 395                                         continue;
 396                         }
 397                         if(fProtonAnalysisBase->GetEtaMode())
 398                         {
 399                                 containerInput[0] =particle->Eta();
 400                         }
 401                         else
 402                         {
 403                                 containerInput[0] = fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz());
 404                         }
 405                         //containerInput[0] = particle->Eta() ;
 406                         containerInput[1] = (Float_t)particle->Pt();
 407                         if (particle->GetPdgCode()==2212)
 408                         {
 409                                 fProtonContainer->Fill(containerInput,kAll,weight);
 410                                 if(ipart<stack->GetNprimary())
 411                                 {
 412                                         fProtonContainer->Fill(containerInput,kPrimary,weight);
 413                                 }
 414                                 else
 415                                 {
 416                                         if(trackflag==1)
 417                                                 fProtonContainer->Fill(containerInput,kFromLambda,weight);
 418                                 }
 419                         }
 420                         if (particle->GetPdgCode()==-2212)
 421                         {
 422                                 fAntiProtonContainer->Fill(containerInput,kAll,weight);
 423                                 if(ipart<stack->GetNprimary())
 424                                 {
 425                                         fAntiProtonContainer->Fill(containerInput,kPrimary,weight);
 426                                 }
 427                                 else
 428                                 {
 429                                         if(trackflag==-1)
 430                                                 fAntiProtonContainer->Fill(containerInput,kFromLambda,weight);
 431                                 }
 432                         }
 433                 }
 434         }
 435 }
 436 //______________________________________________________________________________________________
 437 Int_t AliProtonFeedDownAnalysis::LambdaIsMother(Int_t number, AliStack *stack)
 438 {
 439          if(number<0)
 440                 return 0;
 441         TParticle *particle  = stack->Particle(number);
 442         Int_t motherPDG=0;
 443         Int_t lmother=-1;
 444         lmother=particle->GetFirstMother();
 445         if (lmother<0)
 446                 return 0;
 447         TParticle *mparticle=stack->Particle(lmother);
 448         motherPDG=mparticle->GetPdgCode();
 449         if(motherPDG==3122)
 450                 return 1;
 451         if(motherPDG==-3122)
 452                 return -1;
 453         return 0;
 454 }
 455 //___________________________________________________________________________________________
 456 Float_t AliProtonFeedDownAnalysis::GetWeightforProton(Int_t number,AliStack *stack)
 457 {
 458          if(number<0)
 459                 return 1.0;
 460         TParticle *particle  = stack->Particle(number);
 461         Int_t lmother=particle->GetFirstMother();
 462         TParticle *mparticle=stack->Particle(lmother);
 463         return  fweightfunction->Eval(mparticle->Pt());
 464 }
 465 //______________________________________________________________________________________________
 466 Float_t AliProtonFeedDownAnalysis::GetWeightforLambda(Float_t pt)
 467 {
 468         return  fweightfunction->Eval(pt);
 469 }