PWGLF/RESONANCES/extra/AliAnalysisTaskLambdaStar.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 #if !defined(__CINT__) || defined(__MAKECINT__)
  16
  17 #include <TH1F.h>
  18 #include <TH2F.h>
  19 #include <TH2D.h>
  20 #include <TH3F.h>
  21 #include <TExMap.h>
  22 #include <TVector3.h>
  23 #include <TAxis.h>
  24 #include <TRandom.h>
  25 #include <AliAnalysisTask.h>
  26 #include <AliAnalysisManager.h>
  27
  28 #include <AliAODEvent.h>
  29 #include <AliAODVertex.h>
  30 #include <AliAODv0.h>
  31 #include <AliAODInputHandler.h>
  32
  33 #include "AliAnalysisTaskLambdaStar.h"
  34 #include <AliCentrality.h>
  35 #include <AliPID.h>
  36 #include <AliPIDResponse.h>
  37 #include <AliExternalTrackParam.h>
  38 #include <AliAODTrack.h>
  39 #include <AliVTrack.h>
  40
  41
  42
  43 ClassImp(AliAnalysisTaskLambdaStar)
  44
  45
  46 #endif
  47 //________________________________________________________________________
  48 AliAnalysisTaskLambdaStar::AliAnalysisTaskLambdaStar()
  49   : AliAnalysisTaskSE(),
  50   fkAbsZvertexCut(10.0),
  51   fkCentCut(80.0),
  52   fkKaonMass(0.49366),
  53   fkProMass(0.9382720),
  54   fPIDResponse(0),
  55   fCirc(0),
  56   fCentMin(0),
  57   fCentMax(0),
  58   fNSigma(0),
  59   fNMix(0),
  60   fPatch(0),
  61   fCentPerPatch(0),
  62   fStrong(0),
  63   fResoBuffer(0),
  64   fAOD(0),
  65   fPrimaryVtx(0),
  66   fOutputList(0),
  67   fOutputPrimaries(0),
  68   fGTI(0),fTrackBuffSize(18000),
  69   fHistGoodEvent(0),
  70   fHistEvent(0),
  71   fHistZVertexCent(0),
  72   fPriHistShare(0),
  73   fHistMassPtPKMin(0),
  74   fHistMassPtPbarKPlus(0),
  75   fHistMassPtPKMinMix(0),
  76   fHistMassPtPbarKPlusMix(0),
  77   fHistMassPtPKPlusLS(0),
  78   fHistMassPtPbarKMinLS(0),
  79   fPriHistTPCsignalPos(0),
  80   fPriHistTPCsignalNeg(0),
  81   fPriHistDCAxyYPtPro(0),
  82   fPriHistDCAxyYPtAPro(0),
  83   fPriHistDCAxyYPtKPlus(0),
  84   fPriHistDCAxyYPtKMinus(0)
  85
  86 {
  87   // Dummy constructor
  88   fPrimaryVtxPosition[0]=0;
  89   fPrimaryVtxPosition[1]=0;
  90   fPrimaryVtxPosition[2]=0;
  91 }
  92 //________________________________________________________________________
  93 AliAnalysisTaskLambdaStar::AliAnalysisTaskLambdaStar(const char *name)
  94   : AliAnalysisTaskSE(name),
  95   fkAbsZvertexCut(10.0),
  96   fkCentCut(80.0),
  97   fkKaonMass(0.493),
  98   fkProMass(0.9382720),
  99   fPIDResponse(0),
 100   fCirc(0),
 101   fCentMin(0),
 102   fCentMax(0),
 103   fNSigma(0),
 104   fNMix(0),
 105   fPatch(0),
 106   fCentPerPatch(0),
 107   fStrong(0),
 108   fResoBuffer(0),
 109   fAOD(0),
 110   fPrimaryVtx(0),
 111   fOutputList(0),
 112   fOutputPrimaries(0),
 113   fGTI(0),
 114   fTrackBuffSize(18000),
 115   fHistGoodEvent(0),
 116   fHistEvent(0),
 117   fHistZVertexCent(0),
 118   fPriHistShare(0),
 119   fHistMassPtPKMin(0),
 120   fHistMassPtPbarKPlus(0),
 121   fHistMassPtPKMinMix(0),
 122   fHistMassPtPbarKPlusMix(0),
 123   fHistMassPtPKPlusLS(0),
 124   fHistMassPtPbarKMinLS(0),
 125   fPriHistTPCsignalPos(0),
 126   fPriHistTPCsignalNeg(0),
 127   fPriHistDCAxyYPtPro(0),
 128   fPriHistDCAxyYPtAPro(0),
 129   fPriHistDCAxyYPtKPlus(0),
 130   fPriHistDCAxyYPtKMinus(0)
 131
 132 {
 133   // Constructor
 134   fPrimaryVtxPosition[0]=0;
 135   fPrimaryVtxPosition[1]=0;
 136   fPrimaryVtxPosition[2]=0;
 137
 138   // Define output slots only here
 139   // Output slot #1 writes into a TList container
 140   DefineOutput(1, TList::Class());
 141   DefineOutput(2, TList::Class());
 142
 143 }
 144 //________________________________________________________________________
 145 AliAnalysisTaskLambdaStar::~AliAnalysisTaskLambdaStar() {
 146   // Destructor, go through the data member and delete them
 147
 148   // fPIDResponse is just a pointer to the pid response task,
 149   // we don't create it so we don't delete it. It comes from
 150   // the AliInputEventHandler
 151
 152   if(fResoBuffer){
 153     delete fResoBuffer;
 154     fResoBuffer=0;
 155   }
 156
 157
 158   // The lists containing the histograms
 159   if (fOutputList){
 160     fOutputList->Delete();
 161     delete fOutputList;
 162     fOutputList=0;
 163   }
 164   if (fOutputPrimaries){
 165     fOutputPrimaries->Delete();
 166     delete fOutputPrimaries;
 167     fOutputPrimaries=0;
 168   }
 169
 170
 171   // Array, note the [] with the delete
 172   if (fGTI)
 173     delete[] fGTI;
 174   fGTI=0;
 175
 176 }
 177 //________________________________________________________________________
 178 void AliAnalysisTaskLambdaStar::UserCreateOutputObjects()
 179 {
 180   // Create histograms and other objects and variables
 181   // Called once
 182
 183    // Get the PID response object
 184   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 185   if(!man){AliError("Couldn't get the analysis manager!");}
 186
 187   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 188   if(!inputHandler){AliError("Couldn't get the input handler!");}
 189   fPIDResponse = inputHandler->GetPIDResponse();
 190   if(!fPIDResponse){AliError("Couldn't get the PID response task!");}
 191
 192   // Create the buffer for event mixing
 193   // Standard values are
 194   //  fkZvertexBins(10),
 195   //  fkCentBins(10),
 196   //  fkMixBuff(5),
 197   //  fkPriTrackLim(500),
 198   //  fkV0Lim(50),
 199
 200   fResoBuffer = new ResoBuffer(10,10,fNMix,1200,fkAbsZvertexCut,fkCentCut);
 201
 202   // In AODs, TPC only tracks don't have the pid information stored.
 203   // Also, the TPC only tracks don't have any resolution in the DCAxy
 204   // to distinguish between primaries and secondaries so we need the
 205   // corresponding global track for every TPC only track. The way to do
 206   // this is to just store the pointer to the global track for every id.
 207   fGTI = new AliAODTrack *[fTrackBuffSize]; // Array of pointers
 208
 209   // Create the output list
 210   fOutputList = new TList();
 211   fOutputList->SetOwner();
 212   fOutputPrimaries = new TList();
 213   fOutputPrimaries->SetOwner();
 214
 215   // Invariant mass binning for lambdas
 216   const Int_t nMinvBins = 300;
 217   const Float_t minvLowEdge=1.4, minvHiEdge=1.7;
 218
 219   // Control hist for event cuts
 220   fHistGoodEvent = new TH1F("h1GoodEvent","No of events passing the cuts.",10,-.5,9.5);
 221   fHistEvent = new TH1F("hEvent", "Accepted Event Vs. Centrality", 100,0.0,100);
 222   fHistZVertexCent = new TH2D("fHistZVertexCent"," Vz  Vs Centrality; V_{Z} {cm} ; Centrality {%}", 60, -15, 15,100,0,100);
 223     // 3d y pt mass
 224   fHistMassPtPKMin = new TH2D ("InvMassPtPKMin","M_{inv}pt Vs mass",nMinvBins,minvLowEdge,minvHiEdge,300,0.,30.);
 225   fHistMassPtPbarKPlus = new TH2D ("InvMassPtPKPlus","M_{inv}pt Vs mass",nMinvBins,minvLowEdge,minvHiEdge ,300, 0., 30.0);
 226   fHistMassPtPKMinMix = new TH2D ("InvMassPtPKMinMix","M_{inv}pt Vs mass",nMinvBins,minvLowEdge,minvHiEdge ,300, 0., 30.);
 227   fHistMassPtPbarKPlusMix = new TH2D ("InvMassPtPKPlusMix","M_{inv}pt Vs mass",nMinvBins,minvLowEdge,minvHiEdge , 300, 0., 30.0);
 228   fHistMassPtPKPlusLS = new TH2D ("InvMassPtPKMinLS","M_{inv}pt Vs mass",nMinvBins,minvLowEdge,minvHiEdge ,300, 0.0, 30.);
 229   fHistMassPtPbarKMinLS = new TH2D ("InvMassPtPKPlusLS","M_{inv}pt Vs mass",nMinvBins,minvLowEdge,minvHiEdge ,300,0.0,30.0);
 230
 231   fOutputList->Add(fHistGoodEvent);
 232   fOutputList->Add(fHistEvent);
 233   fOutputList->Add(fHistZVertexCent);
 234   fOutputList->Add(fHistMassPtPKMin);
 235   fOutputList->Add(fHistMassPtPbarKPlus);
 236   fOutputList->Add(fHistMassPtPKMinMix);
 237   fOutputList->Add(fHistMassPtPbarKPlusMix);
 238   fOutputList->Add(fHistMassPtPKPlusLS);
 239   fOutputList->Add(fHistMassPtPbarKMinLS);
 240
 241   // Shared clusters
 242   fPriHistShare = new TH1F ("h1PriShare","Shared clusters, primaries;#shared clusters;counts", 160,0,160);
 243   // dEdx analysis
 244   fPriHistTPCsignalPos = new TH2F ("h2TPCsignalPos","TPC signal for positives;p_{tot};dEdx",400,0,4,1000,0,400);
 245   fPriHistTPCsignalNeg = new TH2F ("h2TPCsignalNeg","TPC signal for negatives;p_{tot};dEdx",400,0.0,4.0,1000,0.0,400.0);
 246
 247  //  Common for all protons - DCA xy distribution to determine primaries, secondaries from weak decay and secondaries from material
 248
 249   fPriHistDCAxyYPtPro = new TH3F ("h3DCAxyYPtPro","DCAxy vs (y,pt) protons",100,-3.,3.,30,-1.5,1.5,100,0.,30);
 250   fPriHistDCAxyYPtAPro = new TH3F ("h3DCAxyYPtAPro","DCAxy vs (y,pt) anti-protons",100,-3.,3.,30,-1.5,1.5,100,0.,30);
 251   fPriHistDCAxyYPtKPlus = new TH3F ("h3DCAxyYPtKPlus","DCAxy vs (y,pt) kplus",100,-3.,3.,30,-1.5,1.5,100,0.,30.);
 252   fPriHistDCAxyYPtKMinus = new TH3F ("h3DCAxyYPtKMinus","DCAxy vs (y,pt) kminus",100,-3.,3.,30,-1.5,1.5,100,0.,30.);
 253   fOutputPrimaries->Add(fPriHistShare);
 254   fOutputPrimaries->Add(fPriHistTPCsignalPos);
 255   fOutputPrimaries->Add(fPriHistTPCsignalNeg);
 256   fOutputPrimaries->Add(fPriHistDCAxyYPtPro);
 257   fOutputPrimaries->Add(fPriHistDCAxyYPtAPro);
 258   fOutputPrimaries->Add(fPriHistDCAxyYPtKPlus);
 259   fOutputPrimaries->Add(fPriHistDCAxyYPtKMinus);
 260
 261   // Post the data
 262   PostData(1, fOutputList);
 263   PostData(2, fOutputPrimaries);
 264
 265
 266 }
 267
 268 //________________________________________________________________________
 269 void AliAnalysisTaskLambdaStar::UserExec(Option_t *)
 270 {
 271   // Main loop
 272   // Called for each event and Fill a control histogram
 273   fHistGoodEvent->Fill(0.0);
 274   // Get the event
 275   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 276   if (!fAOD) {
 277     printf("ERROR: fAOD not available\n");
 278     return;
 279   }
 280
 281   // Fill a control histogram
 282   fHistGoodEvent->Fill(1.0);
 283
 284   // Get the centrality selection
 285   AliCentrality *centrality=NULL;
 286   centrality = fAOD->GetCentrality();
 287   if (!centrality) {
 288     printf ("ERROR: couldn't get the AliCentrality\n");
 289     return;
 290   }
 291
 292   // Fill a control histogram
 293   fHistGoodEvent->Fill(2.0);
 294
 295
 296  if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected())) return;
 297
 298   // Analyze only events using multiplicity in V0 detector (standard)
 299   Float_t centralityPercentile = centrality->GetCentralityPercentileUnchecked("V0M");
 300   if(!centrality->GetCentralityPercentileUnchecked("V0M"))
 301     {
 302     return ;
 303     }
 304
 305   fHistGoodEvent->Fill(3.0);
 306
 307    if ( centralityPercentile <= fCentMin || centralityPercentile > fCentMax){
 308   return;
 309   }
 310
 311    if ( centralityPercentile >= 5 && centralityPercentile <= 20){
 312      ApplyCentralityPatchPbPb2011(centrality);  }
 313
 314      fHistEvent->Fill(centralityPercentile);
 315
 316   // Fill a control histogram
 317      fHistGoodEvent->Fill(4.0);
 318
 319      // Primary vertex, GetPrimaryVertex() returns the "best" reconstructed vertex
 320   fPrimaryVtx = fAOD->GetPrimaryVertex();
 321   if (!fPrimaryVtx){
 322     printf ("ERROR: no primary vertex\n");
 323     return;
 324   }
 325
 326   // Fill a control histogram
 327   fHistGoodEvent->Fill(5.0);
 328
 329   fPrimaryVtx->GetXYZ(fPrimaryVtxPosition);
 330
 331   if (TMath::Abs(fPrimaryVtxPosition[2]) > fkAbsZvertexCut)
 332    return;
 333
 334   // Fill a control histogram
 335   fHistGoodEvent->Fill(6.0);
 336
 337   //Fill  Z-vertex histo
 338   fHistZVertexCent->Fill(fPrimaryVtxPosition[2], centralityPercentile);
 339
 340   // Multiplicity
 341   if (!(fAOD->GetNumberOfTracks())) {
 342     return;
 343   }
 344
 345   // Fill a control histogram
 346   fHistGoodEvent->Fill(7.0);
 347
 348   // Set up the event buffer to store this event
 349   fResoBuffer->ShiftAndAdd(fAOD);
 350
 351   // Reset the reference array to the global tracks..
 352   ResetGlobalTrackReference();
 353
 354   //  AliAODTrack *track=NULL;
 355
 356   //   AliAODTrack* t = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(iTrack));
 357
 358   for (Int_t iTrack=0;iTrack<fAOD->GetNumberOfTracks();iTrack++){
 359
 360     AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
 361     //track = fAOD->GetTrack(iTrack);
 362     if (!track) continue;
 363     // Store the reference of the global tracks
 364     StoreGlobalTrackReference(track);
 365
 366   }
 367   for (Int_t iTrack=0;iTrack<fAOD->GetNumberOfTracks();iTrack++){
 368     AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
 369     // track = fAOD->GetTrack(iTrack);
 370     if (!track) continue;
 371     if(!track->TestFilterBit(1024))
 372       continue;
 373     if(!AcceptTrack(track))
 374       continue;
 375
 376     // Reject tracks with shared clusters
 377     if(!GoodTPCFitMapSharedMap(track))
 378       continue;
 379
 380     // Visualization of TPC dE/dx
 381       FillDedxHist(track);
 382
 383     // Depending on momentum choose pid method
 384     if (track->P() < 0.65){
 385       ProcessTPC(track,fNSigma,fStrong);
 386     }
 387      else if (track->P() < 1.2){
 388        ProcessHybridPro(track,fCirc,fNSigma,fStrong);
 389      }
 390     else if (track->P() < 100.0) {
 391       ProcessHybrid(track,fCirc,fNSigma,fStrong);
 392     }
 393
 394   } // End of loop over primary tracks
 395
 396
 397   // Process real events
 398    ProcessReal( );
 399    ProcessLikeSignBkg();
 400
 401   // Process mixed events
 402    ProcessMixed();
 403
 404   // Post output data.
 405   PostData(1, fOutputList);
 406   PostData(2, fOutputPrimaries);
 407
 408
 409 }
 410
 411
 412 //________________________________________________________________________
 413 void AliAnalysisTaskLambdaStar::ProcessTPC(AliAODTrack* track, Double_t nsig, Bool_t strong){
 414
 415
 416   Double_t nsigmapion = 999, nsigmakaon=999,nsigmaproton=999,nsigmaelectron=999;
 417
 418   nsigmaelectron = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron)) ;
 419   nsigmakaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) ;
 420   nsigmapion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) ;
 421   nsigmaproton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) ;
 422
 423   if( nsigmaelectron  < 3.0  &&  nsigmapion  > 3.0  && nsigmakaon  > 3.0  && nsigmaproton > 3.0 )
 424   return;
 425
 426   if(strong){
 427     if(( nsigmakaon==nsigmapion ) && ( nsigmakaon==nsigmaproton )) return ;
 428     }
 429   //  cout<<"NSigma on TPC Loop = "<<nsig;
 430   if( nsigmakaon <= nsig ){
 431     if (track->Charge()>0){
 432
 433       // Cut .1 cm on DCAxy and fill a histogram
 434       if(goodDCAKaon(track)){
 435
 436         // Add to the  event
 437         fResoBuffer->GetEvt(0)->AddKPlus(track);
 438       }
 439     }
 440     else{
 441       // Cut .1 cm on DCAxy and fill a histogram
 442
 443       if(goodDCAKaon(track)){
 444         // Add to the  event
 445         fResoBuffer->GetEvt(0)->AddKMin(track);
 446       }
 447     }
 448   }
 449
 450    if(strong){
 451    if( ( nsigmaproton==nsigmapion ) && ( nsigmaproton==nsigmakaon )) return ;
 452    }
 453      if( nsigmaproton   <= nsig ){
 454        if (track->Charge()>0){
 455
 456       // Cut .1 cm on DCAxy and fill a histogram
 457          if(goodDCA(track)){
 458         // Add to the  event
 459         fResoBuffer->GetEvt(0)->AddPro(track);
 460       }
 461    }
 462     else{
 463       // Cut .1 cm on DCAxy and fill a histogram
 464       if(goodDCA(track)){
 465         // Add to the  event
 466         fResoBuffer->GetEvt(0)->AddAPro(track);
 467       }
 468     }
 469   }
 470
 471 } // End of void ProcessTPC
 472
 473
 474
 475
 476 //________________________________________________________________________
 477 void AliAnalysisTaskLambdaStar::ProcessHybridPro(AliAODTrack *track, Bool_t circ, Double_t nsig, Bool_t strong){
 478
 479   Double_t nsigmapion = 999, nsigmakaon=999,nsigmaproton=999,nsigmaelectron=999;
 480
 481   nsigmaelectron = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron)) ;
 482   nsigmapion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) ;
 483   nsigmakaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) ;
 484   nsigmaproton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) ;
 485
 486   Double_t nsigmaprotonTOF=999.,nsigmakaonTOF=999.,nsigmapionTOF=999.;
 487   Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
 488
 489   Bool_t fHasTOFPID;
 490
 491   if( nsigmaelectron  < 3.0  &&  nsigmapion  > 3.0  && nsigmakaon  > 3.0  && nsigmaproton > 3.0 )  return;
 492
 493   if(track->GetStatus() & AliVTrack::kTOFpid){
 494     fHasTOFPID=kTRUE;
 495     }
 496     else{
 497       fHasTOFPID=kFALSE;
 498     }
 499
 500
 501   if (fHasTOFPID){
 502     nsigmapionTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion)) ;
 503     nsigmakaonTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon)) ;
 504     nsigmaprotonTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) ;
 505
 506     if(circ){
 507     Double_t d2Proton=nsigmaproton * nsigmaproton + nsigmaprotonTOF * nsigmaprotonTOF;
 508     Double_t d2Kaon=nsigmakaon * nsigmakaon + nsigmakaonTOF * nsigmakaonTOF;
 509     Double_t d2Pion=nsigmapion * nsigmapion + nsigmapionTOF * nsigmapionTOF;
 510
 511     nsigmaTPCTOFkProton  =  TMath::Sqrt(d2Proton);
 512     nsigmaTPCTOFkKaon    =  TMath::Sqrt(d2Kaon);
 513     nsigmaTPCTOFkPion    =  TMath::Sqrt(d2Pion);
 514     }
 515     else{
 516
 517       nsigmaTPCTOFkProton  =  nsigmaprotonTOF;
 518       nsigmaTPCTOFkKaon    =  nsigmakaonTOF;
 519       nsigmaTPCTOFkPion    =  nsigmapionTOF;
 520     }
 521   }
 522   else
 523      {
 524        nsigmaTPCTOFkProton = nsigmaproton;
 525        nsigmaTPCTOFkKaon   = nsigmakaon;
 526        nsigmaTPCTOFkPion   = nsigmapion;
 527      }
 528
 529   if(strong)
 530        {
 531          if(circ){
 532          if( (nsigmaTPCTOFkKaon >= nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon >= nsigmaTPCTOFkProton )) return ;
 533          }
 534          else{
 535          if( (nsigmakaon >=nsigmapion ) && ( nsigmakaon >=nsigmaproton )) return ;
 536          if(fHasTOFPID)
 537            {
 538          if( (nsigmakaonTOF >=nsigmapionTOF ) && ( nsigmakaonTOF >=nsigmaprotonTOF )) return ;
 539            }
 540
 541          }
 542     }
 543
 544   if(!circ) {
 545
 546   if( ( nsigmaTPCTOFkKaon   <= nsig ) && ( nsigmakaon <= nsig ) ){
 547
 548    //
 549
 550   // Distinguish between charges
 551
 552   if (track->Charge() > 0) {
 553
 554         // Cut .1 cm on DCAxy and fill a histogram
 555         if(goodDCAKaon(track)){
 556           // Add to the  event
 557           fResoBuffer->GetEvt(0)->AddKPlus(track);
 558         }
 559       }
 560     else {
 561         // Cut .1 cm on DCAxy and fill a histogram
 562       if(goodDCAKaon(track)){
 563           // add to the  event
 564           fResoBuffer->GetEvt(0)->AddKMin(track);
 565         }
 566     }
 567   }
 568
 569     }
 570
 571   else
 572     {
 573
 574        if( nsigmaTPCTOFkKaon  <= nsig ){
 575
 576          if (track->Charge() > 0) {
 577
 578         // Cut .1 cm on DCAxy and fill a histogram
 579         if(goodDCAKaon(track)){
 580           // Add to the  event
 581           fResoBuffer->GetEvt(0)->AddKPlus(track);
 582         }
 583       }
 584     else {
 585         // Cut .1 cm on DCAxy and fill a histogram
 586       if(goodDCAKaon(track)){
 587           // add to the  event
 588           fResoBuffer->GetEvt(0)->AddKMin(track);
 589         }
 590     }
 591        }
 592     }
 593
 594   if(strong){
 595
 596     if( ( nsigmaproton==nsigmapion ) && ( nsigmaproton==nsigmakaon )) return ;
 597   }
 598
 599   if( nsigmaproton   <= nsig ){
 600
 601   // Distinguish between charges
 602       if (track->Charge() > 0) {
 603
 604         // Cut .1 cm on DCAxy and fill a histogram
 605         if(goodDCA(track)){
 606           // Add to the  event
 607           fResoBuffer->GetEvt(0)->AddPro(track);
 608         }
 609       }
 610     else {
 611         // Cut .1 cm on DCAxy and fill a histogram
 612
 613       if(goodDCA(track)){
 614
 615           // add to the  event
 616
 617           fResoBuffer->GetEvt(0)->AddAPro(track);
 618       }
 619     }
 620
 621   }
 622
 623
 624
 625
 626 } // End of ProcessHybrid
 627 //________________________________________________________________________
 628
 629 void AliAnalysisTaskLambdaStar::ProcessHybrid(AliAODTrack *track, Bool_t circ,Double_t nsig,Bool_t strong){
 630
 631   Double_t nsigmapion = 999, nsigmakaon=999,nsigmaproton=999,nsigmaelectron=999;
 632
 633   nsigmaelectron = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron)) ;
 634   nsigmapion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) ;
 635   nsigmakaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) ;
 636   nsigmaproton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) ;
 637
 638     Double_t nsigmaprotonTOF=999.,nsigmakaonTOF=999.,nsigmapionTOF=999.;
 639     Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
 640
 641     Bool_t fHasTOFPID;
 642
 643     if(track->GetStatus() & AliVTrack::kTOFpid){
 644
 645      fHasTOFPID=kTRUE;
 646     }
 647     else{
 648
 649       fHasTOFPID=kFALSE;
 650     }
 651
 652
 653   if (fHasTOFPID){
 654
 655          nsigmakaonTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon)) ;
 656          nsigmaprotonTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) ;
 657          nsigmapionTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion)) ;
 658
 659          if(circ){
 660
 661            Double_t d2Proton=nsigmaproton * nsigmaproton + nsigmaprotonTOF * nsigmaprotonTOF;
 662            Double_t d2Kaon=nsigmakaon * nsigmakaon + nsigmakaonTOF * nsigmakaonTOF;
 663            Double_t d2Pion=nsigmapion * nsigmapion + nsigmapionTOF * nsigmapionTOF;
 664
 665            nsigmaTPCTOFkProton  =  TMath::Sqrt(d2Proton);
 666            nsigmaTPCTOFkKaon    =  TMath::Sqrt(d2Kaon);
 667            nsigmaTPCTOFkPion    =  TMath::Sqrt(d2Pion);
 668          }
 669     else{
 670
 671       nsigmaTPCTOFkProton  =  nsigmaprotonTOF;
 672       nsigmaTPCTOFkKaon    =  nsigmakaonTOF;
 673       nsigmaTPCTOFkPion    =  nsigmapionTOF;
 674     }
 675   }
 676   else {
 677
 678     nsigmaTPCTOFkProton = nsigmaproton;
 679     nsigmaTPCTOFkKaon   = nsigmakaon;
 680     nsigmaTPCTOFkPion   = nsigmapion;
 681   }
 682
 683   if(strong)
 684        {
 685          if(circ){
 686            if( (nsigmaTPCTOFkKaon >= nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon >= nsigmaTPCTOFkProton )) return ;
 687          }
 688          else{
 689          if( (nsigmakaon >=nsigmapion ) && ( nsigmakaon >=nsigmaproton )) return ;
 690          if(fHasTOFPID){
 691          if( (nsigmakaonTOF >=nsigmapionTOF ) && ( nsigmakaonTOF >=nsigmaprotonTOF )) return ;
 692          }
 693          }
 694     }
 695
 696   if(!circ){
 697
 698       if( ( nsigmaTPCTOFkKaon   <= nsig ) && ( nsigmakaon <= nsig ) ){
 699
 700         if (track->Charge() > 0) {
 701
 702         // Cut .1 cm on DCAxy and fill a histogram
 703         if(goodDCAKaon(track)){
 704           // Add to the  event
 705           fResoBuffer->GetEvt(0)->AddKPlus(track);
 706         }
 707       }
 708     else {
 709         // Cut .1 cm on DCAxy and fill a histogram
 710       if(goodDCAKaon(track)){
 711           // add to the  event
 712           fResoBuffer->GetEvt(0)->AddKMin(track);
 713         }
 714     }
 715   }
 716
 717     }
 718
 719   else
 720     {
 721        if( nsigmaTPCTOFkKaon  <= nsig ){
 722          if (track->Charge() > 0) {
 723            if(goodDCAKaon(track)){
 724            fResoBuffer->GetEvt(0)->AddKPlus(track);
 725         }
 726       }
 727     else {
 728         // Cut .1 cm on DCAxy and fill a histogram
 729       if(goodDCAKaon(track)){
 730           // add to the  event
 731           fResoBuffer->GetEvt(0)->AddKMin(track);
 732         }
 733     }
 734        }
 735     }
 736
 737
 738   //  proton selection
 739
 740   if(strong){
 741          if(circ){
 742            if( (nsigmaTPCTOFkProton >= nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkProton >= nsigmaTPCTOFkKaon )) return ;
 743          }
 744          else{
 745          if( (nsigmaproton >=nsigmapion ) && ( nsigmaproton >=nsigmakaon )) return ;
 746          if(fHasTOFPID){
 747          if( (nsigmaprotonTOF >=nsigmapionTOF ) && ( nsigmaprotonTOF >=nsigmakaonTOF )) return ;
 748          }
 749          }
 750     }
 751
 752   if(!circ){
 753
 754       if( ( nsigmaTPCTOFkProton   <= nsig ) && ( nsigmaproton <= nsig ) ){
 755
 756         if (track->Charge() > 0) {
 757         if(goodDCA(track)){
 758         fResoBuffer->GetEvt(0)->AddPro(track);
 759         }
 760       }
 761     else {
 762
 763       if(goodDCA(track)){
 764       fResoBuffer->GetEvt(0)->AddAPro(track);
 765         }
 766     }
 767   }
 768
 769     }
 770
 771   else
 772     {
 773        if( nsigmaTPCTOFkProton  <= nsig ){
 774          if (track->Charge() > 0) {
 775            if(goodDCA(track)){
 776              fResoBuffer->GetEvt(0)->AddPro(track);
 777         }
 778       }
 779     else {
 780
 781       if(goodDCA(track)){
 782       fResoBuffer->GetEvt(0)->AddAPro(track);
 783
 784       }
 785     }
 786        }
 787     }
 788
 789
 790 } // End of ProcessHybrid
 791
 792
 793   Double_t AliAnalysisTaskLambdaStar::ApplyCentralityPatchPbPb2011( AliCentrality *central){
 794     //This part rejects randomly events such that the centrality gets flat for LHC11h Pb-Pb data
 795     //for 0-5% and 10-20% centrality bin
 796
 797     Double_t cent = (Float_t)(central->GetCentralityPercentile("V0M"));
 798     Double_t rnd_hc, testf, ff, N1, N2;
 799
 800     // cout<<"Centrality patch value"<<fCentPerPatch<<endl;
 801
 802     if(fCentPerPatch==510){
 803
 804       N1 = 1.9404e+06;
 805       N2 = 1.56435e+06;
 806       ff = 5.04167e+06 - 1.49885e+06*cent + 2.35998e+05*cent*cent -1.22873e+04*cent*cent*cent;
 807     }
 808
 809     if(fCentPerPatch==1020){
 810         N1 = 1.56435e+06;
 811         N2 = 4.20e+05;
 812         ff = 1.68062e+08 - 5.19673e+07*cent + 6.4068e+06*cent*cent + 6.4068e+06*cent*cent*cent - 392687*cent*cent*cent*cent - 145.07*cent*cent*cent*cent*cent;
 813     }
 814
 815     testf = ( N2 + (N1-ff) ) / N1;
 816     rnd_hc = gRandom->Rndm();
 817
 818
 819     if (rnd_hc < 0 || rnd_hc > 1 )
 820     {
 821         AliWarning("Wrong Random number generated");
 822         return -999.0;
 823     }
 824
 825     if (rnd_hc < testf)
 826         return cent;
 827     else
 828         return -999.0;
 829 }
 830
 831
 832
 833 //________________________________________________________________________
 834 void AliAnalysisTaskLambdaStar::ProcessReal() {
 835   // Process real events
 836
 837   Int_t iPro,iKminus,iAPro,iKplus;
 838
 839   // Proton K- loop
 840   Int_t nproton = fResoBuffer->GetEvt(0)->GetNPro();
 841   Int_t nkmin = fResoBuffer->GetEvt(0)->GetNKMin();
 842
 843     for (iPro = 0; iPro < nproton; iPro++){
 844     // Skip if unUseIt() entry
 845     if (!fResoBuffer->GetEvt(0)->fProTracks[iPro].UseIt())
 846       continue;
 847     // Kminus loop
 848     for (iKminus=0;iKminus < nkmin;iKminus++){
 849
 850       // Skip if unUseIt() entry
 851       if (!fResoBuffer->GetEvt(0)->fKMinTracks[iKminus].UseIt())
 852         continue;
 853
 854
 855       Double_t  pairrap =  Rapidity(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
 856       Double_t  invmass = MInv(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
 857       Double_t  pairpt  = Pt(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
 858       //      Double_t  ctheta1  = Costheta(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
 859       // Double_t  ctheta2  = Costheta1(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
 860       //Double_t  openang =  OpeningAngle(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
 861
 862       if(TMath::Abs(pairrap) > 0.5) continue;
 863
 864       //if(openang < 0.4) continue;
 865        //  if(TMath::Abs(ctheta1) > 0.8 && TMath::Abs(ctheta2 > 0.8)) continue;
 866        //      cout<<"*****openang"<<openang<<"********"<<ctheta1<<"****"<<ctheta2<<endl;
 867       //cout<<"rap"<<pairrap<<"invmass"<<invmass<<"pairpt"<<pairpt<<endl;
 868
 869       fHistMassPtPKMin->Fill(invmass,pairpt);
 870
 871     }// Kaon loop
 872
 873   }// Proton loop
 874
 875         Int_t npbar   = fResoBuffer->GetEvt(0)->GetNAPro();
 876         Int_t nkplus  = fResoBuffer->GetEvt(0)->GetNKPlus();
 877
 878
 879         for (iAPro = 0; iAPro<npbar; iAPro++){
 880
 881           // Skip if unUseIt() entry
 882
 883           if (!fResoBuffer->GetEvt(0)->fAProTracks[iAPro].UseIt())  continue;
 884
 885           // Kplus loop
 886
 887           for (iKplus=0;iKplus< nkplus;iKplus++){
 888
 889
 890             if (!fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus].UseIt()) continue;
 891
 892             Double_t  pairrap =  Rapidity(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
 893             Double_t  invmass = MInv(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
 894             Double_t  pairpt  = Pt(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
 895             //   Double_t  openang  = OpeningAngle(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
 896             //Double_t  ctheta1  = Costheta(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
 897             //Double_t  ctheta2  = Costheta1(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
 898
 899
 900             if(TMath::Abs(pairrap) > 0.5) continue;
 901
 902
 903             //    if(openang < 0.4) continue;
 904             //   if(TMath::Abs(ctheta1) > 0.8 && TMath::Abs(ctheta2 > 0.8)) continue;
 905             //      cout<<"*****openang"<<openang<<"********"<<ctheta1<<"****"<<ctheta2<<endl;
 906             //cout<<" rap "<<pairrap<<"invmass     "<<invmass<<"    pairpt    "<<pairpt<<endl;
 907
 908             fHistMassPtPbarKPlus->Fill(invmass,pairpt);
 909             // Fill the ThnSparse
 910
 911
 912           }// Kplus loop
 913
 914         }//A Proton loop
 915
 916
 917 }
 918
 919 //________________________________________________________________________
 920 void AliAnalysisTaskLambdaStar::ProcessMixed() {
 921   // Process mixed events
 922
 923   Int_t iPro, iKminus, iAPro, iKplus;
 924
 925   // Int_t nmixed = fResoBuffer->GetMixBuffSize();
 926
 927   //  cout<<"nmixed*******"<<nmixed<<endl;
 928
 929   // Loop over the event buffer
 930   for (UChar_t iMix = 1;iMix<fResoBuffer->GetMixBuffSize();iMix++){
 931
 932     Int_t nproton = fResoBuffer->GetEvt(0)->GetNPro();
 933     Int_t nkmin = fResoBuffer->GetEvt(iMix)->GetNKMin();
 934
 935     for (iPro = 0; iPro < nproton; iPro++){
 936
 937       // Skip if unUseIt() entry
 938       if (!fResoBuffer->GetEvt(0)->fProTracks[iPro].UseIt())
 939         continue;
 940
 941       // Proton loop
 942     for (iKminus=0;iKminus< nkmin;iKminus++){
 943
 944         // Skip if unUseIt() entry
 945         if (!(fResoBuffer->GetEvt(iMix))->fKMinTracks[iKminus].UseIt())
 946           continue;
 947
 948
 949         Double_t  pairrap =  Rapidity(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(iMix)->fKMinTracks[iKminus]);
 950         Double_t  invmass = MInv(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(iMix)->fKMinTracks[iKminus]);
 951         Double_t  pairpt  = Pt(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(iMix)->fKMinTracks[iKminus]);
 952         //      Double_t  openang  = OpeningAngle(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(iMix)->fKMinTracks[iKminus]);
 953
 954         //cout<<invmass<<"****"<<pairpt<<"**********mixed"<<endl;
 955
 956         if(TMath::Abs(pairrap) > 0.5) continue;
 957
 958
 959         //if(openang < 0.4) continue;
 960
 961         fHistMassPtPKMinMix->Fill(invmass,pairpt);
 962
 963    }// Kmin loop
 964
 965     }// Proton loop
 966
 967
 968     Int_t npbar   = fResoBuffer->GetEvt(0)->GetNAPro();
 969
 970     Int_t nkplus  = fResoBuffer->GetEvt(iMix)->GetNKPlus();
 971
 972   //    for (iAPro = 0; iAPro < fResoBuffer->GetEvt(0)->GetNAPro(); iAPro++){
 973
 974     for (iAPro = 0; iAPro < npbar; iAPro++){
 975     // Skip if unUseIt() entry
 976     if (!fResoBuffer->GetEvt(0)->fAProTracks[iAPro].UseIt())
 977       continue;
 978     // Kplus loop
 979     for (iKplus=0;iKplus< nkplus ;iKplus++){
 980
 981       // Skip if unUseIt() entry
 982       if (!fResoBuffer->GetEvt(iMix)->fKPlusTracks[iKplus].UseIt())
 983         continue;
 984
 985        Double_t  pairrap =  Rapidity(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(iMix)->fKPlusTracks[iKplus]);
 986        Double_t  invmass = MInv(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(iMix)->fKPlusTracks[iKplus]);
 987        Double_t  pairpt  = Pt(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(iMix)->fKPlusTracks[iKplus]);
 988
 989        //       Double_t  openang  = OpeningAngle(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(iMix)->fKPlusTracks[iKplus]);
 990
 991
 992         if(TMath::Abs(pairrap) > 0.5) continue;
 993
 994         //if(openang < 0.4) continue;
 995
 996        //cout<<invmass<<"****"<<pairpt<<"**********mixed"<<endl;
 997         fHistMassPtPbarKPlusMix->Fill(invmass,pairpt);
 998
 999     }// AProton loop
1000     }// kplus loop
1001
1002   }// Event buffer loop
1003
1004 }// End of void ProcessMixed
1005 //________________________________________________________________________
1006 void AliAnalysisTaskLambdaStar::ProcessLikeSignBkg() {
1007
1008   Int_t iPro, iKminus, iAPro, iKplus ;
1009
1010   // Proton K+ loop
1011   Int_t nproton = fResoBuffer->GetEvt(0)->GetNPro();
1012   Int_t nkplus = fResoBuffer->GetEvt(0)->GetNKPlus();
1013
1014
1015   for (iPro = 0; iPro < nproton; iPro++){
1016     // Skip if unUseIt() entry
1017     if (!fResoBuffer->GetEvt(0)->fProTracks[iPro].UseIt())
1018       continue;
1019     // Kplus loop
1020     for (iKplus=0;iKplus< nkplus;iKplus++){
1021
1022       // Skip if unUseIt() entry
1023       if (!fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus].UseIt())
1024         continue;
1025
1026       Double_t  pairrap =  Rapidity(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
1027       Double_t  invmass = MInv(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
1028       Double_t  pairpt  = Pt(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
1029       //  Double_t  openang  = OpeningAngle(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
1030       //Double_t  ctheta1  = Costheta(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
1031       //Double_t  ctheta2  = Costheta1(fResoBuffer->GetEvt(0)->fProTracks[iPro], fResoBuffer->GetEvt(0)->fKPlusTracks[iKplus]);
1032
1033       if(TMath::Abs(pairrap) > 0.5) continue;
1034
1035       // if(openang < 0.4) continue;
1036       //cout<<"rap"<<pairrap<<"invmass"<<invmass<<"pairpt"<<pairpt<<endl;
1037
1038       fHistMassPtPKPlusLS->Fill(invmass,pairpt);
1039
1040     }// Kaon loop
1041
1042   }// Proton loop
1043
1044   Int_t npbar = fResoBuffer->GetEvt(0)->GetNAPro();
1045   Int_t nkmin = fResoBuffer->GetEvt(0)->GetNKMin();
1046
1047   for (iAPro = 0; iAPro < npbar; iAPro++){
1048     // Skip if unUseIt() entry
1049     if (!fResoBuffer->GetEvt(0)->fAProTracks[iAPro].UseIt())
1050       continue;
1051     // Kminus loop
1052     for (iKminus=0;iKminus<nkmin;iKminus++){
1053
1054       // Skip if unUseIt() entry
1055       if (!fResoBuffer->GetEvt(0)->fKMinTracks[iKminus].UseIt())
1056         continue;
1057
1058        Double_t  pairrap =  Rapidity(fResoBuffer->GetEvt(0)->fAProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
1059        Double_t  invmass = MInv(fResoBuffer->GetEvt(0)->fAProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
1060        Double_t  pairpt  = Pt(fResoBuffer->GetEvt(0)->fAProTracks[iPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
1061        //  Double_t  openang  = OpeningAngle(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
1062        //Double_t  ctheta1  = Costheta(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
1063        //Double_t  ctheta2  = Costheta1(fResoBuffer->GetEvt(0)->fAProTracks[iAPro], fResoBuffer->GetEvt(0)->fKMinTracks[iKminus]);
1064
1065
1066
1067        if(TMath::Abs(pairrap) > 0.5) continue;
1068
1069        //       if(openang < 0.4) continue;
1070        //cout<<"rap"<<pairrap<<"invmass"<<invmass<<"pairpt"<<pairpt<<endl;
1071
1072        fHistMassPtPbarKMinLS->Fill(invmass,pairpt);
1073
1074        // Fill the ThnSparse
1075
1076
1077     }// Kaon loop
1078
1079   }// Proton loop
1080
1081
1082
1083
1084 } //
1085
1086
1087
1088 Float_t AliAnalysisTaskLambdaStar::MInv(ResoBufferTrack track1, ResoBufferTrack track2){
1089
1090
1091   Float_t e1 = TMath::Sqrt(track1.fP[0]*track1.fP[0] + track1.fP[1]*track1.fP[1] + track1.fP[2]*track1.fP[2] + fkProMass*fkProMass);
1092   Float_t e2 = TMath::Sqrt(track2.fP[0]*track2.fP[0] + track2.fP[1]*track2.fP[1] + track2.fP[2]*track2.fP[2] + fkKaonMass*fkKaonMass);
1093
1094   Double_t massinv =  TMath::Sqrt(((e1+e2) * (e1+e2)) - (((track1.fP[0]+track2.fP[0])*(track1.fP[0]+track2.fP[0])) + ((track1.fP[1]+track2.fP[1])*(track1.fP[1]+track2.fP[1])) + ((track1.fP[2]+track2.fP[2])*(track1.fP[2]+track2.fP[2]))));
1095
1096  return massinv;
1097 }
1098
1099 Float_t AliAnalysisTaskLambdaStar::Costheta(ResoBufferTrack track1, ResoBufferTrack track2){
1100
1101
1102   Double_t massvtx = 1.520 ;
1103   Double_t massp[2];
1104   massp[0] = fkProMass;
1105   massp[1] = fkKaonMass;
1106
1107   Double_t pStar = TMath::Sqrt((massvtx*massvtx-(massp[0]*massp[0])-(massp[1]*massp[1]))*(massvtx*massvtx-(massp[0]*massp[0])-(massp[1]*massp[1]))- (4.* massp[0]*massp[0]*massp[1]*massp[1]))/(2.*massvtx);
1108
1109   Double_t ener =  TMath::Sqrt((massvtx * massvtx) + (((track1.fP[0]+track2.fP[0])*(track1.fP[0]+track2.fP[0])) + ((track1.fP[1]+track2.fP[1])*(track1.fP[1]+track2.fP[1])) + ((track1.fP[2]+track2.fP[2])*(track1.fP[2]+track2.fP[2]))));
1110
1111
1112   Double_t momlam = TMath::Sqrt(((track1.fP[0]+track2.fP[0])*(track1.fP[0]+track2.fP[0])) + ((track1.fP[1]+track2.fP[1])*(track1.fP[1]+track2.fP[1])) + ((track1.fP[2]+track2.fP[2])*(track1.fP[2]+track2.fP[2])));
1113
1114   //  Double_t e=E(pdgvtx);
1115
1116   Double_t beta = momlam/ener;
1117   Double_t gamma = ener/massvtx;
1118   TVector3 mom(track1.fP[0],track1.fP[1],track1.fP[2]);
1119   TVector3 momTot((track1.fP[0] + track2.fP[0]), (track1.fP[1] + track2.fP[1]), (track1.fP[2] + track2.fP[2]));
1120   Double_t QlProng =  mom.Dot(momTot)/momTot.Mag();
1121   Double_t cts = (QlProng/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1122
1123    return  TMath::Cos(cts);
1124
1125 }
1126
1127 Float_t AliAnalysisTaskLambdaStar::Costheta1(ResoBufferTrack track1, ResoBufferTrack track2){
1128
1129
1130   Double_t massvtx = 1.520 ;
1131   Double_t massp[2];
1132   massp[0] = fkProMass;
1133   massp[1] = fkKaonMass;
1134
1135   Double_t pStar = TMath::Sqrt((massvtx*massvtx-(massp[0]*massp[0])-(massp[1]*massp[1]))*(massvtx*massvtx-(massp[0]*massp[0])-(massp[1]*massp[1]))- (4.* massp[0]*massp[0]*massp[1]*massp[1]))/(2.*massvtx);
1136
1137   Double_t ener =  TMath::Sqrt((massvtx * massvtx) + (((track1.fP[0]+track2.fP[0])*(track1.fP[0]+track2.fP[0])) + ((track1.fP[1]+track2.fP[1])*(track1.fP[1]+track2.fP[1])) + ((track1.fP[2]+track2.fP[2])*(track1.fP[2]+track2.fP[2]))));
1138
1139
1140   Double_t momlam = TMath::Sqrt(((track1.fP[0]+track2.fP[0])*(track1.fP[0]+track2.fP[0])) + ((track1.fP[1]+track2.fP[1])*(track1.fP[1]+track2.fP[1])) + ((track1.fP[2]+track2.fP[2])*(track1.fP[2]+track2.fP[2])));
1141
1142   //  Double_t e=E(pdgvtx);
1143
1144   Double_t beta = momlam/ener;
1145   Double_t gamma = ener/massvtx;
1146   TVector3 mom(track2.fP[0],track2.fP[1],track2.fP[2]);
1147   TVector3 momTot((track1.fP[0] + track2.fP[0]), (track1.fP[1] + track2.fP[1]), (track1.fP[2] + track2.fP[2]));
1148   Double_t QlProng =  mom.Dot(momTot)/momTot.Mag();
1149   Double_t cts = (QlProng/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
1150
1151    return  TMath::Cos(cts);
1152
1153 }
1154
1155
1156
1157
1158 //________________________________________________________________________
1159 Float_t AliAnalysisTaskLambdaStar::Rapidity(ResoBufferTrack track1, ResoBufferTrack track2){
1160
1161   Float_t e1 = TMath::Sqrt((track1.fP[0]*track1.fP[0]) + (track1.fP[1]*track1.fP[1]) + (track1.fP[2]*track1.fP[2]) + (fkProMass*fkProMass));
1162   Float_t e2 = TMath::Sqrt((track2.fP[0]*track2.fP[0]) + (track2.fP[1]*track2.fP[1]) + (track2.fP[2]*track2.fP[2]) + (fkKaonMass*fkKaonMass));
1163   Double_t e = e1 + e2;
1164   Double_t pz = track1.fP[2]+ track2.fP[2];
1165
1166   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1167     return 0.5*TMath::Log((e+pz)/(e-pz));
1168   } else { // energy was equal to pz
1169     return -999.;
1170   }
1171
1172 }
1173
1174
1175 Double_t AliAnalysisTaskLambdaStar::OpeningAngle(ResoBufferTrack track1, ResoBufferTrack track2){
1176
1177
1178   Double_t e1e2 = (track1.fP[0]*track2.fP[0]) + (track1.fP[1]*track2.fP[1]) + (track1.fP[2]*track2.fP[2]);
1179   Double_t e2 = TMath::Sqrt((track2.fP[0]*track2.fP[0]) + (track2.fP[1]*track2.fP[1]) + (track2.fP[2]*track2.fP[2]));
1180   Double_t e1 = TMath::Sqrt((track1.fP[0]*track1.fP[0]) + (track1.fP[1]*track1.fP[1]) + (track1.fP[2]*track1.fP[2]));
1181
1182   return TMath::Cos(e1e2/(e1*e2));
1183   //  return   57.296*(TMath::ACos(e1e2/e1*e2));
1184
1185 }
1186
1187
1188
1189 //________________________________________________________________________
1190 Bool_t AliAnalysisTaskLambdaStar::goodDCA(AliAODTrack *track) {
1191
1192   // Get the DCAxy and DCAz. There also exists a TPC only
1193   // impact parameter, but this has not enough resolution
1194   // to discriminate between primaries, secondaries and material
1195
1196   Float_t xy=0.,rap=RapidityProton(track),pt=track->Pt();
1197
1198       xy = DCAxy(track, fAOD);
1199
1200   // Fill the DCAxy histograms
1201   if (track->Charge() > 0){
1202     fPriHistDCAxyYPtPro->Fill(xy,rap,pt);
1203   }
1204   else{
1205     fPriHistDCAxyYPtAPro->Fill(xy,rap,pt);
1206   }
1207   // Do a cut. 0.1 cm shows highest significance for primaries
1208   if (xy>0.1)
1209     return kFALSE;
1210   return kTRUE;
1211 }
1212
1213
1214 Bool_t AliAnalysisTaskLambdaStar::goodDCAKaon(AliAODTrack *track) {
1215
1216   // Get the DCAxy and DCAz. There also exists a TPC only
1217   // impact parameter, but this has not enough resolution
1218   // to discriminate between primaries, secondaries and material
1219
1220   Float_t xy=0.,rap=RapidityKaon(track),pt=track->Pt();
1221
1222         xy = DCAxy(track, fAOD);
1223
1224
1225   // Fill the DCAxy histograms
1226   if (track->Charge() > 0){
1227     fPriHistDCAxyYPtKPlus->Fill(xy,rap,pt);
1228   }
1229   else{
1230     fPriHistDCAxyYPtKMinus->Fill(xy,rap,pt);
1231   }
1232   // Do a cut. 0.1 cm shows highest significance for primaries
1233   if (xy>0.1)
1234     return kFALSE;
1235   return kTRUE;
1236 }
1237 //_______________________________________________________________
1238 Float_t AliAnalysisTaskLambdaStar::RapidityProton(AliAODTrack *track){
1239   // Can't find how to set the assumed mass for the AliAODTrack.
1240   // Same stuff as in AliAODTrack::Y() just with proton mass
1241   Double_t e = TMath::Sqrt(track->P()*track->P() + fkProMass*fkProMass);
1242   Double_t pz = track->Pz();
1243   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1244     return 0.5*TMath::Log((e+pz)/(e-pz));
1245   } else { // energy was equal to pz
1246     return -999.;
1247   }
1248 }
1249
1250
1251 Float_t AliAnalysisTaskLambdaStar::RapidityKaon(AliAODTrack *track){
1252   // Can't find how to set the assumed mass for the AliAODTrack.
1253   // Same stuff as in AliAODTrack::Y() just with proton mass
1254   Double_t e = TMath::Sqrt(track->P()*track->P() + fkKaonMass*fkKaonMass);
1255   Double_t pz = track->Pz();
1256   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1257     return 0.5*TMath::Log((e+pz)/(e-pz));
1258   } else { // energy was equal to pz
1259     return -999.;
1260   }
1261 }
1262
1263
1264 //________________________________________________________________________
1265
1266 //________________________________________________________________________
1267 Float_t AliAnalysisTaskLambdaStar::DCAxy(const AliAODTrack *track, const AliVEvent *evt){
1268   // Note that AliAODTrack::PropagateToDCA() changes the track.
1269   // Don't know whether this is what one wants?
1270   if(!track){
1271     printf("Pointer to track is zero!\n");
1272     return -9999.;
1273   }
1274
1275   // Create an external parameter from the AODtrack
1276   AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1277   // Propagation through the beam pipe would need a correction
1278   // for material, I guess.
1279   if(etp.GetX()>3.) {
1280     printf("This method can be used only for propagation inside the beam pipe\n");
1281     printf("  id: %d, filtermap: %d\n",track->GetID(),track->GetFilterMap());
1282     return -9999.;
1283   }
1284   // Do the propagation
1285   Double_t dca[2]={-9999.,-9999.},covar[3]={0.,0.,0.};
1286   if(!etp.PropagateToDCA(evt->GetPrimaryVertex(),evt->GetMagneticField(),10.,dca,covar)) return -9999.;
1287   // return the DCAxy
1288   return dca[0];
1289 }
1290 //________________________________________________________________________
1291 void AliAnalysisTaskLambdaStar::FillDedxHist(const AliVTrack *track){
1292   // This is for visualization. Fill the the dE/dx histograms
1293   // for all tracks, not only for those, where only the TPC
1294   // is used for PID. Thus avoiding the sharp cut off at a
1295   // momentum of 0.75 GeV/c.
1296
1297   // Positive tracks
1298   if (track->Charge() > 0){
1299
1300     fPriHistTPCsignalPos->Fill(track->GetTPCmomentum(),track->GetTPCsignal());
1301   }
1302
1303   // Negative tracks
1304   else{
1305
1306     fPriHistTPCsignalNeg->Fill(track->GetTPCmomentum(),track->GetTPCsignal());
1307
1308   }
1309 }
1310 //________________________________________________________________________
1311 void AliAnalysisTaskLambdaStar::StoreGlobalTrackReference(AliAODTrack *track){
1312
1313   // Check that the id is positive
1314   if(track->GetID()<0){
1315     //    printf("Warning: track has negative ID: %d\n",track->GetID());
1316     return;
1317   }
1318
1319   // Check id is not too big for buffer
1320   if(track->GetID()>=fTrackBuffSize){
1321     printf("Warning: track ID too big for buffer: ID: %d, buffer %d\n"
1322            ,track->GetID(),fTrackBuffSize);
1323     return;
1324   }
1325
1326   // Warn if we overwrite a track
1327   if(fGTI[track->GetID()]){
1328     // Seems like there are FilterMap 0 tracks
1329     // that have zero TPCNcls, don't store these!
1330     if( (!track->GetFilterMap()) &&
1331         (!track->GetTPCNcls())   )
1332       return;
1333
1334     // Imagine the other way around,
1335     // the zero map zero clusters track
1336     // is stored and the good one wants
1337     // to be added. We ommit the warning
1338     // and just overwrite the 'bad' track
1339     if( fGTI[track->GetID()]->GetFilterMap() ||
1340         fGTI[track->GetID()]->GetTPCNcls()   ){
1341       // If we come here, there's a problem
1342       printf("Warning! global track info already there!");
1343       printf("         TPCNcls track1 %u track2 %u",
1344              (fGTI[track->GetID()])->GetTPCNcls(),track->GetTPCNcls());
1345       printf("         FilterMap track1 %u track2 %u\n",
1346              (fGTI[track->GetID()])->GetFilterMap(),track->GetFilterMap());
1347     }
1348   } // Two tracks same id
1349
1350
1351   // Assign the pointer
1352   (fGTI[track->GetID()]) = track;
1353 }
1354 //________________________________________________________________________
1355 void AliAnalysisTaskLambdaStar::ResetGlobalTrackReference(){
1356   // Sets all the pointers to zero. To be called at
1357   // the beginning or end of an event
1358   for(UShort_t i=0;i<fTrackBuffSize;i++){
1359     fGTI[i]=0;
1360   }
1361 }
1362 //________________________________________________________________________
1363 Bool_t AliAnalysisTaskLambdaStar::AcceptTrack(const AliAODTrack *track){
1364   // Apply additional track cuts
1365
1366
1367   Float_t nCrossed = track->GetTPCClusterInfo(2, 1);
1368   if(nCrossed<70)
1369     return kFALSE;
1370   if(!track->GetTPCNclsF())
1371     return kFALSE; // Note that the AliESDtrackCuts would here return kTRUE
1372   if((nCrossed/track->GetTPCNclsF()) < .8)
1373     return kFALSE;
1374
1375   if(TMath::Abs(track->Eta()) > 0.8) return kFALSE;
1376   if(TMath::Abs(track->Pt()) < 0.2) return kFALSE;
1377
1378   return kTRUE;
1379
1380 }
1381 //________________________________________________________________________
1382
1383 //________________________________________________________________________
1384 Bool_t AliAnalysisTaskLambdaStar::GoodTPCFitMapSharedMap(const AliAODTrack *track){
1385   // Rejects tracks with shared clusters after filling a control histogram
1386   // This overload is used for primaries
1387
1388   // Get the shared maps
1389   const TBits sharedMap = track->GetTPCSharedMap();
1390   // Fill a control histogram
1391   fPriHistShare->Fill(sharedMap.CountBits());
1392   // Reject shared clusters
1393   if((sharedMap.CountBits()) >= 1){
1394     // Bad track, has too many shared clusters!
1395     return kFALSE;
1396   }
1397   return kTRUE;
1398 }
1399 //________________________________________________________________________
1400
1401 //________________________________________________________________________
1402 Float_t AliAnalysisTaskLambdaStar::Pt(ResoBufferTrack track1, ResoBufferTrack track2){
1403
1404   return  TMath::Sqrt((track1.fP[0] + track2.fP[0])*(track1.fP[0] + track2.fP[0])  + (track1.fP[1] + track2.fP[1])*(track1.fP[1] + track2.fP[1]));
1405
1406 }
1407
1408
1409 //________________________________________________________________________
1410 AliAnalysisTaskLambdaStar& AliAnalysisTaskLambdaStar::operator=(const AliAnalysisTaskLambdaStar& atpl)
1411 {
1412   if(this!=&atpl){
1413   // One operation with the atpl to get rid of the warning unused parameter
1414   fPrimaryVtxPosition[0]=atpl.fPrimaryVtxPosition[0];
1415   printf("Assignment operator not implemented\n");
1416   }
1417   return *this;
1418 }
1419 //________________________________________________________________________
1420 void AliAnalysisTaskLambdaStar::Terminate(Option_t *)
1421 {
1422   // Draw result to the screen
1423   // Called once at the end of the query
1424 }
1425 //________________________________________________________________________
1426 //
1427 //
1428 //     Classes in the class AliAnalysisTaskLambdaStar
1429 //         ResoBuffer, ResoBufferEvent, ResoBufferV0 and ResoBufferTrack
1430 //
1431 //________________________________________________________________________
1432 //
1433 //                        ResoBufferTrack
1434 //________________________________________________________________________
1435 AliAnalysisTaskLambdaStar::ResoBufferTrack::ResoBufferTrack():
1436   fID(65535)
1437
1438 {
1439   // Standard constructor, initialize everything with values indicating
1440   // a track that should not be used
1441
1442   // No idea how to initialize the arrays nicely like the fID(65535)..
1443   for (UChar_t i=0;i<3;i++){
1444     fP[i]=-9999.;
1445
1446     for (UChar_t j=0;j<9;j++){
1447       //      fXglobal[j][i]=-9999.;
1448       fXshifted[j][i]=-9999.;
1449     }
1450   }
1451 }
1452 //________________________________________________________________________
1453 AliAnalysisTaskLambdaStar::ResoBufferTrack::ResoBufferTrack(const AliAODTrack *track,const Float_t bfield,const Float_t priVtx[3]):
1454   fID(65535)
1455 {
1456
1457   // Use the function to have the code in one place
1458   Set(track,bfield,priVtx);
1459 }
1460 //________________________________________________________________________
1461
1462 //________________________________________________________________________
1463 void AliAnalysisTaskLambdaStar::ResoBufferTrack::GetShiftedPositionAtShiftedRadii(const AliAODTrack *track, const Float_t bfield, const Float_t priVtx[3]){
1464   // Gets the global position of the track at nine different radii in the TPC
1465   // track is the track you want to propagate
1466   // bfield is the magnetic field of your event
1467   // globalPositionsAtRadii is the array of global positions in the radii and xyz
1468
1469   // Initialize the array to something indicating there was no propagation
1470
1471   for(Int_t i=0;i<9;i++){
1472     for(Int_t j=0;j<3;j++){
1473       fXshifted[i][j]=-9999.;
1474     }
1475   }
1476
1477    // Make a copy of the track to not change parameters of the track
1478   AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1479   //  printf("\nAfter CopyFromVTrack\n");
1480   //  etp.Print();
1481
1482   // The global position of the the track
1483   Double_t xyz[3]={-9999.,-9999.,-9999.};
1484
1485   // Counter for which radius we want
1486   Int_t iR=0;
1487   // The radii at which we get the global positions
1488   // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1489   // Compare squared radii for faster code
1490   Float_t RSquaredWanted[9]={85.*85.,105.*105.,125.*125.,145.*145.,165.*165.,
1491                              185.*185.,205.*205.,225.*225.,245.*245.};
1492   // The shifted radius we are at, squared. Compare squared radii for faster code
1493   Float_t shiftedRadiusSquared=0;
1494
1495   // Propagation is done in local x of the track
1496
1497   for (Float_t x = 58.;x<247.;x+=1.){
1498
1499   // Starts at 83 / Sqrt(2) and goes outwards. 85/Sqrt(2) is the smallest local x
1500     // for global radius 85 cm. x = 245 is the outer radial limit of the TPC when
1501     // the track is straight, i.e. has inifinite pt and doesn't get bent.
1502     // If the track's momentum is smaller than infinite, it will develop a y-component,
1503     // which adds to the global radius
1504
1505     // Stop if the propagation was not succesful. This can happen for low pt tracks
1506     // that don't reach outer radii
1507
1508     if(!etp.PropagateTo(x,bfield))break;
1509
1510     etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1511
1512     // Without shifting the primary vertex to (0.,0.,0.) the next line would just be
1513     // WRONG: globalRadiusSquared = xyz[0]*xyz[0]+xyz[1]*xyz[1];
1514     // but as we shift the primary vertex we want to compare positions at shifted radii.
1515     // I can't draw in ASCII but please take a piece of paper and just visualize it once.
1516
1517     // Changing plus to minus on July10th2012
1518
1519     shiftedRadiusSquared = (xyz[0]-priVtx[0])*(xyz[0]-priVtx[0])
1520                          + (xyz[1]-priVtx[1])*(xyz[1]-priVtx[1]);
1521
1522     // Roughly reached the radius we want
1523     if(shiftedRadiusSquared > RSquaredWanted[iR]){
1524
1525       // Bigger loop has bad precision, we're nearly one centimeter too far,
1526       // go back in small steps.
1527       while (shiftedRadiusSquared>RSquaredWanted[iR]){
1528         x-=.1;
1529         //      printf("propagating to x %5.2f\n",x);
1530         if(!etp.PropagateTo(x,bfield))break;
1531         etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1532         // Added the shifting also here on July11th2012
1533         shiftedRadiusSquared = (xyz[0]-priVtx[0])*(xyz[0]-priVtx[0])
1534                              + (xyz[1]-priVtx[1])*(xyz[1]-priVtx[1]);
1535       }
1536       //      printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",TMath::Sqrt(globalRadiusSquared),x,xyz[0],xyz[1],xyz[2]);
1537       fXshifted[iR][0]=xyz[0]-priVtx[0];
1538       fXshifted[iR][1]=xyz[1]-priVtx[1];
1539       fXshifted[iR][2]=xyz[2]-priVtx[2];
1540       // Indicate we want the next radius
1541       iR+=1;
1542     }
1543     if(iR>=8){
1544       // TPC edge reached
1545       return;
1546     }
1547   }
1548 }
1549 //________________________________________________________________________
1550 void AliAnalysisTaskLambdaStar::ResoBufferTrack::Set(const AliAODTrack *track,const Float_t bfield,const Double_t priVtx[3]){
1551   // Overloaded function
1552   Float_t priVtxPos[3]={(Float_t)priVtx[0],(Float_t)priVtx[1],(Float_t)priVtx[2]};
1553   Set(track,bfield,priVtxPos);
1554 }
1555
1556
1557 void AliAnalysisTaskLambdaStar::ResoBufferTrack::SetP(const AliAODTrack *track){
1558   fP[0] = track->Px();
1559   fP[1] = track->Py();
1560   fP[2] = track->Pz();
1561
1562 }
1563 //________________________________________________________________________
1564 void AliAnalysisTaskLambdaStar::ResoBufferTrack::Set(const AliAODTrack *track,const Float_t bfield,const Float_t priVtx[3]){
1565
1566   // Set the ID, a good ID also indicates to use the track
1567   if(track->GetID() >=0){
1568     // global tracks, i.e. v0 daughters
1569     fID = track->GetID();
1570   }
1571   else {
1572     // e.g. tpc only tracks, i.e. primary protons
1573     fID = -track->GetID()-1;
1574
1575   }
1576   // Set the momentum
1577
1578    track->PxPyPz(fP);
1579
1580
1581   GetShiftedPositionAtShiftedRadii(track,bfield,priVtx);
1582
1583 }
1584 //________________________________________________________________________
1585 AliAnalysisTaskLambdaStar::ResoBufferTrack::ResoBufferTrack(const ResoBufferTrack& fbt):
1586   fID(fbt.fID)
1587  {
1588   // Copy constructor
1589
1590   for (UChar_t i=0;i<3;i++){
1591     fP[i]=fbt.fP[i];
1592     for (UChar_t j=0;j<9;j++){
1593       //      fXglobal[j][i]=fbt.fXglobal[j][i];
1594       fXshifted[j][i]=fbt.fXshifted[j][i];
1595     }
1596   }
1597  }
1598 //________________________________________________________________________
1599 AliAnalysisTaskLambdaStar::ResoBufferTrack& AliAnalysisTaskLambdaStar::ResoBufferTrack::operator=(const ResoBufferTrack& fbt){
1600   // Assignment operator, from wikipedia :)
1601
1602   // Protect against self-assignment
1603   if(this != &fbt){
1604     fID = fbt.fID;
1605     for (UChar_t i=0;i<3;i++){
1606       fP[i]=fbt.fP[i];
1607       for (UChar_t j=0;j<9;j++){
1608         //      fXglobal[j][i]=fbt.fXglobal[j][i];
1609         fXshifted[j][i]=fbt.fXshifted[j][i];
1610       }
1611     }
1612   }
1613   // By convention, always return *this (Could it be the convention is called
1614   return *this;
1615 }
1616 //________________________________________________________________________
1617
1618
1619
1620
1621
1622
1623 //                        ResoBufferEvent
1624 //________________________________________________________________________
1625 AliAnalysisTaskLambdaStar::ResoBufferEvent::ResoBufferEvent():
1626   fPriTrackLim(0)
1627   ,fProTracks(0),fAProTracks(0)
1628   ,fKPlusTracks(0),fKMinTracks(0)
1629   ,fNProTracks(0),fNAProTracks(0),fNKPlusTracks(0),fNKMinTracks(0)
1630   ,fBfield(-9999.)
1631 {
1632   // Standard constructor, all pointer to zero
1633   fPriVtxPos[0]=-9999.;
1634   fPriVtxPos[1]=-9999.;
1635   fPriVtxPos[2]=-9999.;
1636
1637   printf("This constructor has zero size in the arrays!\n");
1638 }
1639 //________________________________________________________________________
1640 AliAnalysisTaskLambdaStar::ResoBufferEvent::ResoBufferEvent(const UShort_t priTrackBuff,const Double_t bfield,const Double_t priVtxPos[3]):
1641   fPriTrackLim(priTrackBuff)
1642   ,fProTracks(new ResoBufferTrack[fPriTrackLim])
1643   ,fAProTracks(new ResoBufferTrack[fPriTrackLim])
1644   ,fKPlusTracks(new ResoBufferTrack[fPriTrackLim])
1645   ,fKMinTracks(new ResoBufferTrack[fPriTrackLim])
1646   ,fNProTracks(0),fNAProTracks(0),fNKPlusTracks(0),fNKMinTracks(0)
1647   ,fBfield(-bfield)
1648
1649 {
1650   // Constructor.
1651   fPriVtxPos[0] = priVtxPos[0]; // This is some old C++
1652   fPriVtxPos[1] = priVtxPos[1];
1653   fPriVtxPos[2] = priVtxPos[2];
1654 }
1655 //________________________________________________________________________
1656 AliAnalysisTaskLambdaStar::ResoBufferEvent::ResoBufferEvent(const UShort_t priTrackBuff):
1657   fPriTrackLim(priTrackBuff)
1658   ,fProTracks(new ResoBufferTrack[fPriTrackLim])
1659   ,fAProTracks(new ResoBufferTrack[fPriTrackLim])
1660   ,fKPlusTracks(new ResoBufferTrack[fPriTrackLim])
1661   ,fKMinTracks(new ResoBufferTrack[fPriTrackLim])
1662   ,fNProTracks(0),fNAProTracks(0),fNKPlusTracks(0),fNKMinTracks(0)
1663   ,fBfield(-9999.)
1664
1665 {
1666   // Constructor. fBfield and fPriVtxPos not needed yet, can be set later.
1667   fPriVtxPos[0] = -9999.; // This is C++03
1668   fPriVtxPos[1] = -9999.;
1669   fPriVtxPos[2] = -9999.;
1670
1671   //  printf("constructed eventwith NBgLam: %u\n",fNBgLamTracks);
1672 }
1673 //________________________________________________________________________
1674 AliAnalysisTaskLambdaStar::ResoBufferEvent::ResoBufferEvent(const ResoBufferEvent &fbe):
1675   fPriTrackLim(fbe.GetPriTrackLim())
1676   ,fProTracks(new ResoBufferTrack[fPriTrackLim])
1677   ,fAProTracks(new ResoBufferTrack[fPriTrackLim])
1678   ,fKPlusTracks(new ResoBufferTrack[fPriTrackLim])
1679   ,fKMinTracks(new ResoBufferTrack[fPriTrackLim])
1680   ,fNProTracks(fbe.GetNPro()),fNAProTracks(fbe.GetNAPro())
1681   ,fNKPlusTracks(fbe.GetNKPlus()),fNKMinTracks(fbe.GetNKMin())
1682   ,fBfield(fbe.GetBfield())
1683 {
1684   // Copy constructor
1685   fbe.GetVtxPos(fPriVtxPos);
1686   // Avoid to much creation and deletion of objects
1687   UShort_t i;
1688   // Copy the primary tracks
1689   for (i=0;i<fPriTrackLim;i++){
1690     fProTracks[i]=fbe.fProTracks[i];
1691     fAProTracks[i]=fbe.fAProTracks[i];
1692     fKPlusTracks[i]=fbe.fKPlusTracks[i];
1693     fKMinTracks[i]=fbe.fKMinTracks[i];
1694   }
1695
1696 }
1697 //________________________________________________________________________
1698 AliAnalysisTaskLambdaStar::ResoBufferEvent& AliAnalysisTaskLambdaStar::ResoBufferEvent::operator=(const ResoBufferEvent &fbe){
1699   // Assignment operator
1700
1701   // Protect against self-assignment
1702   if(this!=&fbe){
1703     // Well, we use arrays of a constant size to avoid
1704     // excessive memory allocation and won't give this up.
1705     // So we'll only copy as much as fits on the left side
1706     // from the right side.
1707     // DON'T COPY THE ARRAY SIZES fV0Lim AND fPriTrackLim !!!
1708     if(fPriTrackLim < fbe.GetPriTrackLim() ){
1709       // AliWarning(Form("Trying to assign too big event (buffer %d/%d) to"
1710       //                    " this (buffer %d/%d). Only partially copying.",
1711       //                    fbe.GetPriTrackLim(),
1712       //                    fPriTrackLim,));
1713       printf("Trying to assign too big event (buffer %d) to"
1714                     " this (buffer %d. Only partially copying.\n",
1715              fbe.GetPriTrackLim(),
1716              fPriTrackLim);
1717     }
1718     // Always start with the easy stuff :)
1719     fbe.GetVtxPos(fPriVtxPos);
1720     fBfield = fbe.GetBfield();
1721     // Number of tracks is minimum of array size of 'this'
1722     // and the number of tracks from the right side
1723     fNProTracks = TMath::Min(fPriTrackLim,fbe.GetNPro());
1724     fNAProTracks = TMath::Min(fPriTrackLim,fbe.GetNAPro());
1725
1726     fNKPlusTracks = TMath::Min(fPriTrackLim,fbe.GetNKPlus());
1727     fNKMinTracks = TMath::Min(fPriTrackLim,fbe.GetNKMin());
1728
1729     // Avoid creation and deletion of 'i' for every loop
1730     UShort_t i;
1731     // Copy primary tracks. No need to set a 'bad track'
1732     // flag for the entries above GetNPro() (...) as
1733     // above everything is bad by definition.
1734     // Protons
1735     for (i=0;i<GetNPro();i++)
1736       fProTracks[i]=fbe.fProTracks[i];
1737     // Anti-protons
1738     for (i=0;i<GetNAPro();i++)
1739       fAProTracks[i]=fbe.fAProTracks[i];
1740
1741     // Kplus
1742     for (i=0;i<GetNPro();i++){
1743       fKPlusTracks[i]=fbe.fKPlusTracks[i];
1744     }
1745     // Anti-lambdas
1746     for (i=0;i<GetNKMin();i++){
1747       fKMinTracks[i]=fbe.fKMinTracks[i];
1748     }
1749
1750   }
1751   return *this;
1752 }
1753 //________________________________________________________________________
1754 AliAnalysisTaskLambdaStar::ResoBufferEvent::~ResoBufferEvent(){
1755   // Destructor
1756
1757   // Delete the arrays of tracks,
1758   // note the [] with the delete
1759   if(fProTracks){
1760     delete[] fProTracks;
1761     fProTracks=0;
1762   }
1763   if(fAProTracks){
1764     delete[] fAProTracks;
1765     fAProTracks=0;
1766   }
1767   if(fKPlusTracks){
1768     delete[] fKPlusTracks;
1769     fKPlusTracks=0;
1770   }
1771   if(fKMinTracks){
1772     delete[] fKMinTracks;
1773     fKMinTracks=0;
1774   }
1775
1776 }
1777
1778 //________________________________________________________________________
1779 void AliAnalysisTaskLambdaStar::ResoBufferEvent::Reset(const Double_t bfield, const Double_t priVtxPos[3]){
1780
1781  // Reset the old event, i.e., make clear 'here is no info'
1782   // by setting the 'number of stored ...' to zero
1783   fNProTracks=0;
1784   fNAProTracks=0;
1785   fNKPlusTracks=0;
1786   fNKMinTracks=0;
1787
1788   // And set the new event properties
1789   fBfield = bfield;
1790   fPriVtxPos[0]=priVtxPos[0];
1791   fPriVtxPos[1]=priVtxPos[1];
1792   fPriVtxPos[2]=priVtxPos[2];
1793 }
1794 //________________________________________________________________________
1795 void AliAnalysisTaskLambdaStar::ResoBufferEvent::AddPro(const AliAODTrack *track){
1796   // Add a proton to this event
1797
1798   // Check whether there is still space in the array
1799   if(fNProTracks > fPriTrackLim-1){
1800     // AliWarning(Form("Cannot add proton, array size (%d) too small"
1801     //              ,fPriTrackLim));
1802     printf("Cannot add proton, array size (%d) too small\n"
1803                     ,fPriTrackLim);
1804     return;
1805   }
1806
1807   fProTracks[fNProTracks].Set(track,fBfield,fPriVtxPos);
1808   //  fProTracks[fNProTracks].SetP(track);
1809
1810   fNProTracks++;
1811   //  printf("Added proton %d/%d\n",fNProTracks,fPriTrackLim);
1812
1813 }
1814 //________________________________________________________________________
1815 void AliAnalysisTaskLambdaStar::ResoBufferEvent::AddAPro(const AliAODTrack *track){
1816   // Add a anti-proton to this event
1817
1818
1819   if(fNAProTracks > fPriTrackLim-1){
1820
1821     printf("Cannot add anti-proton, array size (%d) too small\n"
1822                     ,fPriTrackLim);
1823     return;
1824   }
1825   // Add the V0 at the end of the array
1826
1827   fAProTracks[fNAProTracks].Set(track,fBfield,fPriVtxPos);
1828
1829   //fAProTracks[fNAProTracks].SetP(track);
1830
1831   fNAProTracks++;
1832 }
1833 //________________________________________________________________________
1834
1835 void AliAnalysisTaskLambdaStar::ResoBufferEvent::AddKPlus(const AliAODTrack *track){
1836   // Add a kplus to this event
1837
1838   // Check whether there is still space in the array
1839   if(fNKPlusTracks > fPriTrackLim-1){
1840
1841     printf("Cannot add kplus, array size (%d) too small\n"
1842                     ,fPriTrackLim);
1843     return;
1844   }
1845
1846   fKPlusTracks[fNKPlusTracks].Set(track,fBfield,fPriVtxPos);
1847   //fKPlusTracks[fNKPlusTracks].SetP(track);
1848   fNKPlusTracks++;
1849
1850
1851 }
1852 //________________________________________________________________________
1853 void AliAnalysisTaskLambdaStar::ResoBufferEvent::AddKMin(const AliAODTrack *track){
1854   // Add a anti-proton to this event
1855
1856
1857   if(fNKMinTracks > fPriTrackLim-1){
1858
1859     printf("Cannot add kminus, array size (%d) too small\n"
1860                     ,fPriTrackLim);
1861     return;
1862   }
1863   // Add the V0 at the end of the array
1864   fKMinTracks[fNKMinTracks].Set(track,fBfield,fPriVtxPos);
1865   //fKMinTracks[fNKMinTracks].SetP(track);
1866   fNKMinTracks++;
1867 }
1868
1869 //________________________________________________________________________
1870
1871
1872 //________________________________________________________________________
1873 //
1874 //                        ResoBuffer
1875 //________________________________________________________________________
1876 AliAnalysisTaskLambdaStar::ResoBuffer::ResoBuffer() :
1877   fkZvertexBins(0),
1878   fkCentBins(0),
1879   fkMixBuffSize(0),
1880   fkPriTrackLim(0),
1881   fZvertexAxis(0),
1882   fCentAxis(0),
1883   fCurEvt(0),
1884   fEC(0)
1885 {
1886   // Dummy constructor, create arrays with zero size
1887   // Note that some data member are constant, you
1888   // won't be able to create the ResoBuffer first with this
1889   // constructor and then set the appropiate size.
1890
1891 }
1892 //________________________________________________________________________
1893 AliAnalysisTaskLambdaStar::ResoBuffer::ResoBuffer(const UChar_t ZvertexBins,const UChar_t CentBins,const UChar_t MixBuff,const UShort_t PriTrackLim, const Float_t AbsZvertexCut,const Float_t CentCut) :
1894   fkZvertexBins(ZvertexBins),
1895   fkCentBins(CentBins),
1896   fkMixBuffSize(MixBuff),
1897   fkPriTrackLim(PriTrackLim),
1898   fZvertexAxis(new TAxis(fkZvertexBins,-AbsZvertexCut,AbsZvertexCut)),
1899   fCentAxis(new TAxis (fkCentBins,0.0,CentCut)),
1900   fCurEvt(new ResoBufferEvent *[fkMixBuffSize]),
1901   fEC(new ResoBufferEvent ***[fkZvertexBins])
1902 {
1903   // Constructor, creates at once all events with all tracks
1904   //  printf ("Creating with pritracklim %d and v0lim %d\n",fkPriTrackLim,fkV0Lim);
1905
1906   // Create the array step by step
1907   // Bins in z of the primary vertex position. Do this as
1908   // the detector looks different from a different z coordinate
1909   for (UChar_t iZBin=0;iZBin<fkZvertexBins;iZBin++){
1910     fEC[iZBin] = new ResoBufferEvent **[fkCentBins];
1911     // Bins in centrality
1912     for (UChar_t iCentBin=0;iCentBin<fkCentBins;iCentBin++){
1913       fEC[iZBin][iCentBin] = new ResoBufferEvent *[fkMixBuffSize];
1914       // The number of events to keep for one mixing class
1915       for(UChar_t iMixBuff=0;iMixBuff<fkMixBuffSize;iMixBuff++){
1916         // Create an event to hold the info for mixing
1917         fEC[iZBin][iCentBin][iMixBuff] = new ResoBufferEvent(fkPriTrackLim);
1918       }
1919     }
1920   }
1921 }
1922 //________________________________________________________________________
1923 AliAnalysisTaskLambdaStar::ResoBuffer::ResoBuffer(const AliAnalysisTaskLambdaStar::ResoBuffer &fb) :
1924   fkZvertexBins(fb.fkZvertexBins),
1925   fkCentBins(fb.fkCentBins),
1926   fkMixBuffSize(fb.fkMixBuffSize),
1927   fkPriTrackLim(fb.fkPriTrackLim),
1928   fZvertexAxis(new TAxis(*(fb.fZvertexAxis))),
1929   fCentAxis(new TAxis (*(fb.fCentAxis))),
1930   fCurEvt(new ResoBufferEvent *[fkMixBuffSize]),
1931   fEC(new ResoBufferEvent ***[fkZvertexBins])
1932 {
1933   // Copy constructor. Linux complains not having this and
1934   // compiling this task with aliroot
1935
1936   printf("ResoBuffer ctor not tested yet, be cautious\n");
1937
1938   // Create the array step by step
1939   // Bins in z of the primary vertex position. Do this as
1940   // the detector looks different from a different z coordinate
1941   for (UChar_t iZBin=0;iZBin<fkZvertexBins;iZBin++){
1942     fEC[iZBin] = new ResoBufferEvent **[fkCentBins];
1943     // Bins in centrality
1944     for (UChar_t iCentBin=0;iCentBin<fkCentBins;iCentBin++){
1945       fEC[iZBin][iCentBin] = new ResoBufferEvent *[fkMixBuffSize];
1946       // The number of events to keep for one mixing class
1947       for(UChar_t iMixBuff=0;iMixBuff<fkMixBuffSize;iMixBuff++){
1948         // Create an event to hold the info for mixing
1949         fEC[iZBin][iCentBin][iMixBuff] = new ResoBufferEvent(*(fb.fEC[iZBin][iCentBin][iMixBuff]));
1950       }
1951     }
1952   }
1953 }
1954 //________________________________________________________________________
1955 AliAnalysisTaskLambdaStar::ResoBuffer& AliAnalysisTaskLambdaStar::ResoBuffer::operator=(const AliAnalysisTaskLambdaStar::ResoBuffer& fb){
1956   //Assignment operator
1957   if(this!=&fb){
1958     printf("ResoBuffer assignment operator not implemented\n");
1959   }
1960   return *this;
1961
1962 }
1963 //________________________________________________________________________
1964 AliAnalysisTaskLambdaStar::ResoBuffer::~ResoBuffer(){
1965   // Destructor
1966   // The axes to fin the correct bins
1967   if(fZvertexAxis){
1968     delete fZvertexAxis;
1969     fZvertexAxis=0;
1970   }
1971   if(fCentAxis){
1972     delete fCentAxis;
1973     fCentAxis=0;
1974   }
1975   // fCurEvt is an array of pointer
1976   if(fCurEvt){
1977     delete[] fCurEvt;
1978     fCurEvt=0;
1979   }
1980   // Delete all the events and the pointer to them
1981   for (UChar_t iZBin=0;iZBin<fkZvertexBins;iZBin++){
1982     for (UChar_t iCentBin=0;iCentBin<fkCentBins;iCentBin++){
1983       for(UChar_t iMixBuff=0;iMixBuff<fkMixBuffSize;iMixBuff++){
1984         if(fEC[iZBin][iCentBin][iMixBuff]){
1985           delete fEC[iZBin][iCentBin][iMixBuff];
1986           fEC[iZBin][iCentBin][iMixBuff]=0;
1987         }
1988       }
1989       if(fEC[iZBin][iCentBin]){
1990         delete fEC[iZBin][iCentBin];
1991         fEC[iZBin][iCentBin]=0;
1992       }
1993     }
1994     if(fEC[iZBin]){
1995       delete fEC[iZBin];
1996       fEC[iZBin]=0;
1997     }
1998   }
1999   if(fEC){
2000     delete fEC;
2001     fEC=0;
2002   }
2003 }
2004 //________________________________________________________________________
2005 void AliAnalysisTaskLambdaStar::ResoBuffer::ShiftAndAdd(AliAODEvent *evt){
2006   // Shift the events in the appropiate centrality / zvertex bin and set the
2007   // current event pointer correctly
2008   Double_t priVtxPos[3];
2009   evt->GetPrimaryVertex()->GetXYZ(priVtxPos);
2010   //  printf("Mag field: %f\n",evt->GetMagneticField());
2011   ShiftAndAdd(evt->GetMagneticField(),
2012               priVtxPos,
2013               evt->GetCentrality()->GetCentralityPercentileUnchecked("V0M"));
2014 }
2015 //________________________________________________________________________
2016 void AliAnalysisTaskLambdaStar::ResoBuffer::ShiftAndAdd(const Double_t bfield,const Double_t priVtxPos[3],const Float_t centrality){
2017
2018   // Shift the events in the appropiate centrality / zvertex bin and set the
2019   // current event pointer correctly
2020
2021   // Find the correct centrality/zvertex bin
2022
2023   const UChar_t ZvertexBin = fZvertexAxis->FindFixBin(priVtxPos[2]) - 1; // -1 for array starting at 0
2024
2025   const UChar_t CentBin = fCentAxis->FindFixBin(centrality) - 1;// -1 for array starting at 0
2026
2027   // The new current event is the old last event
2028   fCurEvt[0] = fEC[ZvertexBin][CentBin][fkMixBuffSize-1];
2029
2030   // Shift the pointer, starting from the back
2031   UChar_t iMix;
2032   for(iMix=fkMixBuffSize-1;iMix>0;iMix--){
2033     fEC[ZvertexBin][CentBin][iMix] = fEC[ZvertexBin][CentBin][iMix-1];
2034   }
2035   // And reset the zero'th one
2036   fEC[ZvertexBin][CentBin][0] = fCurEvt[0];
2037   fEC[ZvertexBin][CentBin][0]->Reset(bfield,priVtxPos);
2038   // Also set the pointer to the other events..
2039   for (iMix=1;iMix<fkMixBuffSize;iMix++){
2040     fCurEvt[iMix] = fEC[ZvertexBin][CentBin][iMix];
2041   }
2042 }