PWGHF/vertexingHF/AliAnalysisTaskSELc2pK0sfromAODtracks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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  * appeuear 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
  16 /* $Id$ */
  17
  18 //
  19 //
  20 //               Lc->pK0s analysis code
  21 //
  22 //  Input: AOD
  23 //  Output: TTree or THnSparse (mass vs pT vs Centrality)
  24 //
  25 //  Cuts:
  26 //  TTree: very loose cut
  27 //  THnSparse: One THnSparse is created per cut. One cut is specified by
  28 //  an array of bits, each bit corresponds to a cut in "Cut" function.
  29 //  Use "AddCutStream" function to add a cut.
  30 //
  31 //-------------------------------------------------------------------------
  32 //
  33 //                 Authors: Y.S Watanabe(a)
  34 //  (a) CNS, the University of Tokyo
  35 //  Contatcs: wyosuke@cns.s.u-tokyo.ac.jp
  36 //-------------------------------------------------------------------------
  37
  38 #include <TSystem.h>
  39 #include <TParticle.h>
  40 #include <TParticlePDG.h>
  41 #include <TH1F.h>
  42 #include <TH1F.h>
  43 #include <TH2F.h>
  44 #include <THnSparse.h>
  45 #include <TLorentzVector.h>
  46 #include <TTree.h>
  47 #include "TROOT.h"
  48 #include <TDatabasePDG.h>
  49 #include <AliAnalysisDataSlot.h>
  50 #include <AliAnalysisDataContainer.h>
  51 #include "AliStack.h"
  52 #include "AliMCEvent.h"
  53 #include "AliAnalysisManager.h"
  54 #include "AliAODMCHeader.h"
  55 #include "AliAODHandler.h"
  56 #include "AliLog.h"
  57 #include "AliExternalTrackParam.h"
  58 #include "AliAODVertex.h"
  59 #include "AliESDVertex.h"
  60 #include "AliAODRecoDecay.h"
  61 #include "AliAODRecoDecayHF.h"
  62 #include "AliAODRecoCascadeHF.h"
  63 #include "AliESDtrack.h"
  64 #include "AliAODTrack.h"
  65 #include "AliAODv0.h"
  66 #include "AliAODcascade.h"
  67 #include "AliAODMCParticle.h"
  68 #include "AliAnalysisTaskSE.h"
  69 #include "AliAnalysisTaskSELc2pK0sfromAODtracks.h"
  70 #include "AliPIDResponse.h"
  71 #include "AliPIDCombined.h"
  72 #include "AliTOFPIDResponse.h"
  73 #include "AliAODPidHF.h"
  74 #include "AliInputEventHandler.h"
  75 #include "AliESDtrackCuts.h"
  76 #include "AliNeutralTrackParam.h"
  77 #include "AliKFParticle.h"
  78 #include "AliKFVertex.h"
  79 #include "AliExternalTrackParam.h"
  80 #include "AliESDtrack.h"
  81 #include "AliCentrality.h"
  82 #include "AliVertexerTracks.h"
  83
  84 using std::cout;
  85 using std::endl;
  86
  87 ClassImp(AliAnalysisTaskSELc2pK0sfromAODtracks)
  88
  89 //__________________________________________________________________________
  90 AliAnalysisTaskSELc2pK0sfromAODtracks::AliAnalysisTaskSELc2pK0sfromAODtracks() :
  91   AliAnalysisTaskSE(),
  92   fUseMCInfo(kFALSE),
  93   fOutput(0),
  94   fOutputAll(0),
  95   fListCuts(0),
  96   fCEvents(0),
  97   fHTrigger(0),
  98   fHCentrality(0),
  99   fProdCuts(0),
 100   fAnalCuts(0),
 101   fIsEventSelected(kFALSE),
 102   fWriteVariableTree(kFALSE),
 103   fVariablesTree(0),
 104   fIspp(kFALSE),
 105   fIspA(kFALSE),
 106   fIsAA(kFALSE),
 107   fIsMB(kFALSE),
 108   fIsSemi(kFALSE),
 109   fIsCent(kFALSE),
 110   fIsINT7(kFALSE),
 111   fIsEMC7(kFALSE),
 112   fCandidateVariables(),
 113   fVtx1(0),
 114   fV1(0),
 115   fBzkG(0),
 116   fCentrality(0),
 117   fTriggerCheck(0),
 118   fHistoLcK0SpMass(),
 119   fHistoBachPt(0),
 120   fHistod0Bach(0),
 121   fHistod0V0(0),
 122   fHistod0d0(0),
 123   fHistoV0CosPA(0),
 124   fHistoProbProton(0),
 125   fHistoDecayLength(0),
 126   fHistoK0SMass(0)
 127 {
 128   //
 129   // Default Constructor.
 130   //
 131 }
 132
 133 //___________________________________________________________________________
 134 AliAnalysisTaskSELc2pK0sfromAODtracks::AliAnalysisTaskSELc2pK0sfromAODtracks(const Char_t* name,
 135                                                                              AliRDHFCutsLctopK0sfromAODtracks* prodCuts,
 136                                                                              AliRDHFCutsLctopK0sfromAODtracks* analCuts,
 137                                                                              Bool_t writeVariableTree) :
 138   AliAnalysisTaskSE(name),
 139   fUseMCInfo(kFALSE),
 140   fOutput(0),
 141   fOutputAll(0),
 142   fListCuts(0),
 143   fCEvents(0),
 144   fHTrigger(0),
 145   fHCentrality(0),
 146   fProdCuts(prodCuts),
 147   fAnalCuts(analCuts),
 148   fIsEventSelected(kFALSE),
 149   fWriteVariableTree(writeVariableTree),
 150   fVariablesTree(0),
 151   fIspp(kFALSE),
 152   fIspA(kFALSE),
 153   fIsAA(kFALSE),
 154   fIsMB(kFALSE),
 155   fIsSemi(kFALSE),
 156   fIsCent(kFALSE),
 157   fIsINT7(kFALSE),
 158   fIsEMC7(kFALSE),
 159   fCandidateVariables(),
 160   fVtx1(0),
 161   fV1(0),
 162   fBzkG(0),
 163   fCentrality(0),
 164   fTriggerCheck(0),
 165   fHistoLcK0SpMass(),
 166   fHistoBachPt(0),
 167   fHistod0Bach(0),
 168   fHistod0V0(0),
 169   fHistod0d0(0),
 170   fHistoV0CosPA(0),
 171   fHistoProbProton(0),
 172   fHistoDecayLength(0),
 173   fHistoK0SMass(0)
 174 {
 175   //
 176   // Constructor. Initialization of Inputs and Outputs
 177   //
 178   Info("AliAnalysisTaskSELc2pK0sfromAODtracks","Calling Constructor");
 179
 180   DefineOutput(1,TList::Class());  //conters
 181   DefineOutput(2,TList::Class());
 182   if(writeVariableTree){
 183     DefineOutput(3,TTree::Class());  //My private output
 184   }else{
 185     DefineOutput(3,TList::Class());  //conters
 186   }
 187 }
 188
 189 //___________________________________________________________________________
 190 AliAnalysisTaskSELc2pK0sfromAODtracks::~AliAnalysisTaskSELc2pK0sfromAODtracks() {
 191   //
 192   // destructor
 193   //
 194   Info("~AliAnalysisTaskSELc2pK0sfromAODtracks","Calling Destructor");
 195
 196   if (fOutput) {
 197     delete fOutput;
 198     fOutput = 0;
 199   }
 200
 201   if (fOutputAll) {
 202     delete fOutputAll;
 203     fOutputAll = 0;
 204   }
 205
 206   if (fListCuts) {
 207     delete fListCuts;
 208     fListCuts = 0;
 209   }
 210
 211
 212   if (fProdCuts) {
 213     delete fProdCuts;
 214     fProdCuts = 0;
 215   }
 216
 217   if (fAnalCuts) {
 218     delete fAnalCuts;
 219     fAnalCuts = 0;
 220   }
 221
 222   if (fVariablesTree) {
 223     delete fVariablesTree;
 224     fVariablesTree = 0;
 225   }
 226
 227 }
 228
 229 //_________________________________________________
 230 void AliAnalysisTaskSELc2pK0sfromAODtracks::Init() {
 231   //
 232   // Initialization
 233   //
 234   //
 235
 236   fIsEventSelected=kFALSE;
 237
 238   if (fDebug > 1) AliInfo("Init");
 239
 240   fListCuts = new TList();
 241   fListCuts->SetOwner();
 242   fListCuts->SetName("ListCuts");
 243   fListCuts->Add(new AliRDHFCutsLctopK0sfromAODtracks(*fProdCuts));
 244   fListCuts->Add(new AliRDHFCutsLctopK0sfromAODtracks(*fAnalCuts));
 245   PostData(2,fListCuts);
 246
 247   return;
 248 }
 249
 250 //_________________________________________________
 251 void AliAnalysisTaskSELc2pK0sfromAODtracks::UserExec(Option_t *)
 252 {
 253   //
 254   // UserExec
 255   //
 256
 257   if (!fInputEvent) {
 258     AliError("NO EVENT FOUND!");
 259     return;
 260   }
 261   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 262
 263   fCEvents->Fill(1);
 264
 265   //------------------------------------------------
 266   // First check if the event has proper vertex and B
 267   //------------------------------------------------
 268
 269   fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 270   if (!fVtx1) return;
 271
 272   Double_t pos[3],cov[6];
 273   fVtx1->GetXYZ(pos);
 274   fVtx1->GetCovarianceMatrix(cov);
 275   fV1 = new AliESDVertex(pos,cov,100.,100,fVtx1->GetName());
 276
 277   fBzkG = (Double_t)aodEvent->GetMagneticField();
 278   AliKFParticle::SetField(fBzkG);
 279   if (TMath::Abs(fBzkG)<0.001) {
 280     delete fV1;
 281     return;
 282   }
 283   fCEvents->Fill(2);
 284
 285   //------------------------------------------------
 286   // Event selection setting
 287   //------------------------------------------------
 288   if (fUseMCInfo)
 289     fAnalCuts->SetTriggerClass("");
 290
 291   if ( !fUseMCInfo && fIspp) {
 292     fAnalCuts->SetUsePhysicsSelection(kTRUE);
 293     fAnalCuts->SetTriggerClass("");
 294     fAnalCuts->SetTriggerMask(AliVEvent::kMB);
 295   }
 296
 297   if ( !fUseMCInfo && fIspA) {
 298     fAnalCuts->SetUsePhysicsSelection(kTRUE);
 299     fAnalCuts->SetTriggerClass("");
 300     fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
 301   }
 302   if ( !fUseMCInfo && fIsAA) {
 303     fAnalCuts->SetUsePhysicsSelection(kTRUE);
 304     fAnalCuts->SetTriggerClass("");
 305     fAnalCuts->SetTriggerMask(AliVEvent::kMB | AliVEvent::kSemiCentral | AliVEvent::kCentral);
 306   }
 307
 308   //------------------------------------------------
 309   // Event selection
 310   //------------------------------------------------
 311   Bool_t fIsTriggerNotOK = fAnalCuts->IsEventRejectedDueToTrigger();
 312   if(!fIsTriggerNotOK) fCEvents->Fill(3);
 313   fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent); // better to initialize before CheckEventSelection call
 314   if(!fIsEventSelected) {
 315     delete fV1;
 316     return;
 317   }
 318   fCEvents->Fill(4);
 319
 320   fIsMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB)==(AliVEvent::kMB);
 321   fIsSemi=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kSemiCentral)==(AliVEvent::kSemiCentral);
 322   fIsCent=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kCentral)==(AliVEvent::kCentral);
 323   fIsINT7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kINT7)==(AliVEvent::kINT7);
 324   fIsEMC7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kEMC7)==(AliVEvent::kEMC7);
 325   fTriggerCheck = fIsMB+2*fIsSemi+4*fIsCent+8*fIsINT7+16*fIsEMC7;
 326   if(fIsMB) fHTrigger->Fill(1);
 327   if(fIsSemi) fHTrigger->Fill(2);
 328   if(fIsCent) fHTrigger->Fill(3);
 329   if(fIsINT7) fHTrigger->Fill(4);
 330   if(fIsEMC7) fHTrigger->Fill(5);
 331   if(fIsMB|fIsSemi|fIsCent) fHTrigger->Fill(7);
 332   if(fIsINT7|fIsEMC7) fHTrigger->Fill(8);
 333   if(fIsMB&fIsSemi) fHTrigger->Fill(10);
 334   if(fIsMB&fIsCent) fHTrigger->Fill(11);
 335   if(fIsINT7&fIsEMC7) fHTrigger->Fill(12);
 336
 337   if ( !fUseMCInfo && (fIsAA || fIspA)) {
 338     AliCentrality *cent = aodEvent->GetCentrality();
 339     fCentrality = cent->GetCentralityPercentile("V0M");
 340   }
 341   fHCentrality->Fill(fCentrality);
 342
 343   //------------------------------------------------
 344   // Check if the event has v0 candidate
 345   //------------------------------------------------
 346   //Int_t nv0 = aodEvent->GetNumberOfV0s();
 347   fCEvents->Fill(5);
 348
 349   //------------------------------------------------
 350   // MC analysis setting
 351   //------------------------------------------------
 352   //  TClonesArray *mcArray = 0;
 353   //  AliAODMCHeader *mcHeader=0;
 354   //  if (fUseMCInfo) {
 355   //    // MC array need for maching
 356   //    mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 357   //    if (!mcArray) {
 358   //      AliError("Could not find Monte-Carlo in AOD");
 359   //      return;
 360   //    }
 361   //    fCEvents->Fill(6); // in case of MC events
 362   //
 363   //    // load MC header
 364   //    mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 365   //    if (!mcHeader) {
 366   //      AliError("AliAnalysisTaskSELc2pK0sfromAODtracks::UserExec: MC header branch not found!\n");
 367   //      return;
 368   //    }
 369   //    fCEvents->Fill(7); // in case of MC events
 370   //
 371   //    Double_t zMCVertex = mcHeader->GetVtxZ();
 372   //    if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()) {
 373   //      AliDebug(2,Form("Event rejected: abs(zVtxMC)=%f > fAnalCuts->GetMaxVtxZ()=%f",zMCVertex,fAnalCuts->GetMaxVtxZ()));
 374   //      return;
 375   //    } else {
 376   //      fCEvents->Fill(17); // in case of MC events
 377   //    }
 378   //    }
 379
 380   //------------------------------------------------
 381   // Main analysis done in this function
 382   //------------------------------------------------
 383   MakeAnalysis(aodEvent);
 384
 385
 386   PostData(1,fOutput);
 387   if(fWriteVariableTree){
 388     PostData(3,fVariablesTree);
 389   }else{
 390     PostData(3,fOutputAll);
 391   }
 392
 393   fIsEventSelected=kFALSE;
 394
 395   delete fV1;
 396   return;
 397 }
 398
 399 //________________________________________ terminate ___________________________
 400 void AliAnalysisTaskSELc2pK0sfromAODtracks::Terminate(Option_t*)
 401 {
 402   // The Terminate() function is the last function to be called during
 403   // a query. It always runs on the client, it can be used to present
 404   // the results graphically or save the results to file.
 405
 406   //AliInfo("Terminate","");
 407   AliAnalysisTaskSE::Terminate();
 408
 409   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 410   if (!fOutput) {
 411     AliError("fOutput not available");
 412     return;
 413   }
 414
 415   if(!fWriteVariableTree){
 416     fOutputAll = dynamic_cast<TList*> (GetOutputData(3));
 417     if (!fOutputAll) {
 418       AliError("fOutputAll not available");
 419       return;
 420     }
 421   }
 422
 423   return;
 424 }
 425
 426 //___________________________________________________________________________
 427 void AliAnalysisTaskSELc2pK0sfromAODtracks::UserCreateOutputObjects()
 428 {
 429   //
 430   // UserCreateOutputObject
 431   //
 432   //AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
 433
 434   //------------------------------------------------
 435   // output object setting
 436   //------------------------------------------------
 437   fOutput = new TList();
 438   fOutput->SetOwner();
 439   fOutput->SetName("chist0");
 440   DefineGeneralHistograms(); // define general histograms
 441   PostData(1,fOutput);
 442
 443   DefineTreeVariables();
 444   if (fWriteVariableTree) {
 445     PostData(3,fVariablesTree);
 446   }else{
 447     fOutputAll = new TList();
 448     fOutputAll->SetOwner();
 449     fOutputAll->SetName("anahisto");
 450     DefineAnalysisHistograms(); // define general histograms
 451     PostData(3,fOutputAll);
 452   }
 453
 454
 455   return;
 456 }
 457
 458 //-------------------------------------------------------------------------------
 459 void AliAnalysisTaskSELc2pK0sfromAODtracks::MakeAnalysis
 460 (
 461  AliAODEvent *aodEvent
 462  )
 463 {
 464   //
 465   // Main Analysis part
 466   //
 467
 468   Int_t nV0s= aodEvent->GetNumberOfV0s();
 469   if (nV0s==0) {
 470     return;
 471   }
 472   Int_t nTracks= aodEvent->GetNumberOfTracks();
 473   if (nTracks==0) {
 474     return;
 475   }
 476
 477   //------------------------------------------------
 478   // V0 loop
 479   //------------------------------------------------
 480   for (Int_t iv0 = 0; iv0<nV0s; iv0++) {
 481     AliAODv0 *v0 = aodEvent->GetV0(iv0);
 482     if(!v0) continue;
 483     if(!fProdCuts->SingleV0Cuts(v0,fVtx1)) continue;
 484
 485     AliAODTrack *cptrack =  (AliAODTrack*)(v0->GetDaughter(0));
 486     AliAODTrack *cntrack =  (AliAODTrack*)(v0->GetDaughter(1));
 487
 488     //------------------------------------------------
 489     // track loop
 490     //------------------------------------------------
 491     for (Int_t itrk = 0; itrk<nTracks; itrk++) {
 492       AliAODTrack *trk = (AliAODTrack*)aodEvent->GetTrack(itrk);
 493       if(trk->GetID()<0) continue;
 494       if(!fProdCuts->SingleTrkCuts(trk)) continue;
 495
 496       Int_t cpid = cptrack->GetID();
 497       Int_t cnid = cntrack->GetID();
 498       Int_t lpid = trk->GetID();
 499       if((cpid==lpid)||(cnid==lpid)) continue;
 500
 501       if(!fProdCuts->SelectWithRoughCuts(v0,trk)) continue;
 502
 503       AliAODVertex *secVert = ReconstructSecondaryVertex(v0,trk,aodEvent);
 504       if(!secVert) continue;
 505
 506       AliAODRecoCascadeHF *lcobj = MakeCascadeHF(v0,trk,aodEvent,secVert);
 507       if(!lcobj) {
 508         continue;
 509       }
 510
 511       FillROOTObjects(lcobj);
 512
 513       lcobj->GetSecondaryVtx()->RemoveDaughters();
 514       lcobj->UnsetOwnPrimaryVtx();
 515       delete lcobj;lcobj=NULL;
 516       delete secVert;
 517     }
 518   }
 519
 520 }
 521
 522 ////-------------------------------------------------------------------------------
 523 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillROOTObjects(AliAODRecoCascadeHF *lcobj)
 524 {
 525   //
 526   // Fill histograms or tree depending on fWriteVariableTree
 527   //
 528
 529   AliAODTrack *trk = lcobj->GetBachelor();
 530   AliAODv0 *v0 = lcobj->Getv0();
 531
 532   fCandidateVariables[ 0] = lcobj->InvMassLctoK0sP();
 533   fCandidateVariables[ 1] = lcobj->Px();
 534   fCandidateVariables[ 2] = lcobj->Py();
 535   fCandidateVariables[ 3] = lcobj->Pz();
 536   fCandidateVariables[ 4] = v0->MassK0Short();
 537   fCandidateVariables[ 5] = lcobj->PxProng(0);
 538   fCandidateVariables[ 6] = lcobj->PyProng(0);
 539   fCandidateVariables[ 7] = lcobj->PzProng(0);
 540   fCandidateVariables[ 8] = lcobj->PxProng(1);
 541   fCandidateVariables[ 9] = lcobj->PyProng(1);
 542   fCandidateVariables[10] = lcobj->PzProng(1);
 543   fCandidateVariables[11] = fVtx1->GetX();
 544   fCandidateVariables[12] = fVtx1->GetY();
 545   fCandidateVariables[13] = fVtx1->GetZ();
 546   fCandidateVariables[14] = fCentrality;
 547   fCandidateVariables[15] = lcobj->DecayLengthXY();
 548
 549   Double_t nSigmaTPCpr=-9999.;
 550   Double_t nSigmaTOFpr=-9999.;
 551   Double_t probProton=-9999.;
 552   if(fAnalCuts->GetIsUsePID())
 553     {
 554       fAnalCuts->GetPidHF()->GetnSigmaTPC(trk,4,nSigmaTPCpr);
 555       fAnalCuts->GetPidHF()->GetnSigmaTOF(trk,4,nSigmaTOFpr);
 556       probProton = fAnalCuts->GetProtonProbabilityTPCTOF(trk);
 557       fCandidateVariables[16] = nSigmaTPCpr;
 558       fCandidateVariables[17] = nSigmaTOFpr;
 559       fCandidateVariables[18] = probProton;
 560     }
 561
 562   if(fWriteVariableTree)
 563     fVariablesTree->Fill();
 564   else{
 565     for(Int_t ic=0;ic < fAnalCuts->GetNCutsArray(); ic++)
 566       {
 567         fAnalCuts->SetCutsFromArray(ic);
 568         if(fAnalCuts->IsSelected(lcobj,AliRDHFCuts::kCandidate))
 569           {
 570             Double_t cont[3];
 571             cont[0] = lcobj->InvMassLctoK0sP();
 572             cont[1] = lcobj->Pt();
 573             cont[2] = fCentrality;
 574             fHistoLcK0SpMass[ic]->Fill(cont);
 575           }
 576       }
 577     fHistoBachPt->Fill(trk->Pt());
 578     fHistod0Bach->Fill(lcobj->Getd0Prong(0));
 579     fHistod0V0->Fill(lcobj->Getd0Prong(1));
 580     fHistod0d0->Fill(lcobj->Getd0Prong(0)*lcobj->Getd0Prong(1));
 581     fHistoV0CosPA->Fill(lcobj->CosV0PointingAngle());
 582     fHistoProbProton->Fill(probProton);
 583     fHistoDecayLength->Fill(lcobj->DecayLengthXY());
 584     fHistoK0SMass->Fill(v0->MassK0Short());
 585   }
 586   return;
 587 }
 588
 589 ////-------------------------------------------------------------------------------
 590 void AliAnalysisTaskSELc2pK0sfromAODtracks::DefineTreeVariables()
 591 {
 592   //
 593   // Define tree variables
 594   //
 595
 596   const char* nameoutput = GetOutputSlot(3)->GetContainer()->GetName();
 597   fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
 598   Int_t nVar = 19;
 599   fCandidateVariables = new Float_t [nVar];
 600   TString * fCandidateVariableNames = new TString[nVar];
 601
 602   fCandidateVariableNames[ 0]="InvMassLc2pK0s";
 603   fCandidateVariableNames[ 1]="LcPx";
 604   fCandidateVariableNames[ 2]="LcPy";
 605   fCandidateVariableNames[ 3]="LcPz";
 606   fCandidateVariableNames[ 4]="massK0Short";
 607   fCandidateVariableNames[ 5]="V0Px";
 608   fCandidateVariableNames[ 6]="V0Py";
 609   fCandidateVariableNames[ 7]="V0Pz";
 610   fCandidateVariableNames[ 8]="BachPx";
 611   fCandidateVariableNames[ 9]="BachPy";
 612   fCandidateVariableNames[10]="BachPz";
 613   fCandidateVariableNames[11]="PrimVertx";
 614   fCandidateVariableNames[12]="PrimVerty";
 615   fCandidateVariableNames[13]="PrimVertz";
 616   fCandidateVariableNames[14]="Centrality";
 617   fCandidateVariableNames[15]="DecayLengthXY";
 618   fCandidateVariableNames[16]="nSigmaTPCpr";
 619   fCandidateVariableNames[17]="nSigmaTOFpr";
 620   fCandidateVariableNames[18]="probProton";
 621
 622   for (Int_t ivar=0; ivar<nVar; ivar++) {
 623     fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
 624   }
 625
 626   return;
 627 }
 628
 629 ////__________________________________________________________________________
 630 void  AliAnalysisTaskSELc2pK0sfromAODtracks::DefineGeneralHistograms() {
 631   //
 632   // This is to define general histograms
 633   //
 634
 635   fCEvents = new TH1F("fCEvents","conter",18,-0.5,17.5);
 636   fCEvents->SetStats(kTRUE);
 637   fCEvents->GetXaxis()->SetBinLabel(1,"X1");
 638   fCEvents->GetXaxis()->SetBinLabel(2,"Analyzed events");
 639   fCEvents->GetXaxis()->SetBinLabel(3,"AliAODVertex exists");
 640   fCEvents->GetXaxis()->SetBinLabel(4,"TriggerOK");
 641   fCEvents->GetXaxis()->SetBinLabel(5,"IsEventSelected");
 642   fCEvents->GetXaxis()->SetBinLabel(6,"CascadesHF exists");
 643   fCEvents->GetXaxis()->SetBinLabel(7,"MCarray exists");
 644   fCEvents->GetXaxis()->SetBinLabel(8,"MCheader exists");
 645   fCEvents->GetXaxis()->SetBinLabel(9,"triggerClass!=CINT1");
 646   fCEvents->GetXaxis()->SetBinLabel(10,"triggerMask!=kAnyINT");
 647   fCEvents->GetXaxis()->SetBinLabel(11,"triggerMask!=kAny");
 648   fCEvents->GetXaxis()->SetBinLabel(12,"vtxTitle.Contains(Z)");
 649   fCEvents->GetXaxis()->SetBinLabel(13,"vtxTitle.Contains(3D)");
 650   fCEvents->GetXaxis()->SetBinLabel(14,"vtxTitle.Doesn'tContain(Z-3D)");
 651   fCEvents->GetXaxis()->SetBinLabel(15,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
 652   fCEvents->GetXaxis()->SetBinLabel(16,"!IsEventSelected");
 653   fCEvents->GetXaxis()->SetBinLabel(17,"triggerMask!=kAnyINT || triggerClass!=CINT1");
 654   fCEvents->GetXaxis()->SetBinLabel(18,Form("zVtxMC<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
 655   //fCEvents->GetXaxis()->SetTitle("");
 656   fCEvents->GetYaxis()->SetTitle("counts");
 657
 658   fHTrigger = new TH1F("fHTrigger","counter",18,-0.5,17.5);
 659   fHTrigger->SetStats(kTRUE);
 660   fHTrigger->GetXaxis()->SetBinLabel(1,"X1");
 661   fHTrigger->GetXaxis()->SetBinLabel(2,"kMB");
 662   fHTrigger->GetXaxis()->SetBinLabel(3,"kSemiCentral");
 663   fHTrigger->GetXaxis()->SetBinLabel(4,"kCentral");
 664   fHTrigger->GetXaxis()->SetBinLabel(5,"kINT7");
 665   fHTrigger->GetXaxis()->SetBinLabel(6,"kEMC7");
 666   //fHTrigger->GetXaxis()->SetBinLabel(7,"Space");
 667   fHTrigger->GetXaxis()->SetBinLabel(8,"kMB|kSemiCentral|kCentral");
 668   fHTrigger->GetXaxis()->SetBinLabel(9,"kINT7|kEMC7");
 669   fHTrigger->GetXaxis()->SetBinLabel(11,"kMB&kSemiCentral");
 670   fHTrigger->GetXaxis()->SetBinLabel(12,"kMB&kCentral");
 671   fHTrigger->GetXaxis()->SetBinLabel(13,"kINT7&kEMC7");
 672
 673   fHCentrality = new TH1F("fHCentrality","conter",100,0.,100.);
 674
 675   fOutput->Add(fCEvents);
 676   fOutput->Add(fHTrigger);
 677   fOutput->Add(fHCentrality);
 678
 679   return;
 680 }
 681 //__________________________________________________________________________
 682 void  AliAnalysisTaskSELc2pK0sfromAODtracks::DefineAnalysisHistograms()
 683 {
 684   //
 685   // Define analyis histograms
 686   //
 687
 688   //------------------------------------------------
 689   // Basic histogram
 690   //------------------------------------------------
 691   Int_t bins_base[3]=           {80                             ,20             ,10};
 692   Double_t xmin_base[3]={2.286-0.2,0            ,0.00};
 693   Double_t xmax_base[3]={2.286+0.2,20.  ,100};
 694
 695   for(Int_t i=0;i<fAnalCuts->GetNCutsArray();i++)
 696     {
 697       fHistoLcK0SpMass.push_back(new THnSparseF(Form("fHistoLcK0SpMass_Cut%d",i),Form("Cut%d",i),3,bins_base,xmin_base,xmax_base));
 698       fOutputAll->Add(fHistoLcK0SpMass[i]);
 699     }
 700
 701   //------------------------------------------------
 702   // checking histograms
 703   //------------------------------------------------
 704   fHistoBachPt = new TH1F("fHistoBachPt","Bachelor p_{T}",100,0.0,5.0);
 705   fOutputAll->Add(fHistoBachPt);
 706   fHistod0Bach = new TH1F("fHistod0Bach","Bachelor d_{0}",100,-0.5,0.5);
 707   fOutputAll->Add(fHistod0Bach);
 708   fHistod0V0 = new TH1F("fHistod0V0","V_{0} d_{0}",100,-0.5,0.5);
 709   fOutputAll->Add(fHistod0V0);
 710   fHistod0d0 = new TH1F("fHistod0d0","d_{0} d_{0}",100,-0.5,0.5);
 711   fOutputAll->Add(fHistod0d0);
 712   fHistoV0CosPA=new TH1F("fHistoV0CosPA","V0->Second vertex cospa",100,-1.,1.0);
 713   fOutputAll->Add(fHistoV0CosPA);
 714   fHistoProbProton=new TH1F("fHistoProbProton","ProbProton",100,0.,1.0);
 715   fOutputAll->Add(fHistoProbProton);
 716   fHistoDecayLength=new TH1F("fHistoDecayLength","Decay Length",100,-0.1,0.1);
 717   fOutputAll->Add(fHistoDecayLength);
 718   fHistoK0SMass=new TH1F("fHistoK0SMass","K0S mass",100,0.497-0.05,0.497+0.05);
 719   fOutputAll->Add(fHistoK0SMass);
 720
 721   return;
 722 }
 723
 724 //________________________________________________________________________
 725 AliAODRecoCascadeHF* AliAnalysisTaskSELc2pK0sfromAODtracks::MakeCascadeHF(AliAODv0 *v0, AliAODTrack *part, AliAODEvent * aod, AliAODVertex *secVert)
 726 {
 727   //
 728   // Create AliAODRecoCascadeHF object from the argument
 729   //
 730
 731   if(!v0) return 0x0;
 732   if(!part) return 0x0;
 733   if(!aod) return 0x0;
 734
 735   //------------------------------------------------
 736   // PrimaryVertex
 737   //------------------------------------------------
 738   AliAODVertex *primVertexAOD;
 739   Bool_t unsetvtx = kFALSE;
 740   if(fIspp){
 741     primVertexAOD = CallPrimaryVertex(v0,part,aod);
 742     if(!primVertexAOD){
 743       primVertexAOD = fVtx1;
 744     }else{
 745       unsetvtx = kTRUE;
 746     }
 747   }else{
 748     primVertexAOD = fVtx1;
 749   }
 750   if(!primVertexAOD) return 0x0;
 751   Double_t posprim[3]; primVertexAOD->GetXYZ(posprim);
 752
 753   //------------------------------------------------
 754   // DCA between tracks
 755   //------------------------------------------------
 756   AliESDtrack *esdtrack = new AliESDtrack((AliVTrack*)part);
 757
 758   AliNeutralTrackParam *trackV0=NULL;
 759   const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
 760   if(trackVV0)  trackV0 = new AliNeutralTrackParam(trackVV0);
 761
 762   Double_t xdummy, ydummy;
 763   Double_t dca = esdtrack->GetDCA(trackV0,fBzkG,xdummy,ydummy);
 764
 765
 766   //------------------------------------------------
 767   // Propagate all tracks to the secondary vertex and calculate momentum there
 768   //------------------------------------------------
 769
 770   Double_t d0z0bach[2],covd0z0bach[3];
 771   if(sqrt(pow(secVert->GetX(),2)+pow(secVert->GetY(),2))<1.){
 772     part->PropagateToDCA(secVert,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
 773     trackV0->PropagateToDCA(secVert,fBzkG,kVeryBig);
 774   }else{
 775     part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
 776     trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig);
 777   }
 778   Double_t momv0_new[3]={-9999,-9999,-9999.};
 779   trackV0->GetPxPyPz(momv0_new);
 780
 781   Double_t px[2],py[2],pz[2];
 782   px[0] = part->Px(); py[0] = part->Py(); pz[0] = part->Pz();
 783   px[1] = momv0_new[0]; py[1] = momv0_new[1]; pz[1] = momv0_new[2];
 784
 785   //------------------------------------------------
 786   // d0
 787   //------------------------------------------------
 788   Double_t d0[3],d0err[3];
 789
 790   part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
 791   d0[0]= d0z0bach[0];
 792   d0err[0] = TMath::Sqrt(covd0z0bach[0]);
 793
 794   Double_t d0z0v0[2],covd0z0v0[3];
 795   trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0v0,covd0z0v0);
 796   d0[1]= d0z0v0[0];
 797   d0err[1] = TMath::Sqrt(covd0z0v0[0]);
 798
 799   //------------------------------------------------
 800   // Create AliAODRecoCascadeHF
 801   //------------------------------------------------
 802   Short_t charge = part->Charge();
 803   AliAODRecoCascadeHF *theCascade = new AliAODRecoCascadeHF(secVert,charge,px,py,pz,d0,d0err,dca);
 804   if(!theCascade)
 805     {
 806       if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
 807       if(esdtrack) delete esdtrack;
 808       if(trackV0) delete trackV0;
 809       return 0x0;
 810     }
 811   theCascade->SetOwnPrimaryVtx(primVertexAOD);
 812   UShort_t id[2]={(UShort_t)part->GetID(),(UShort_t)trackV0->GetID()};
 813   theCascade->SetProngIDs(2,id);
 814
 815   theCascade->GetSecondaryVtx()->AddDaughter(part);
 816   theCascade->GetSecondaryVtx()->AddDaughter(v0);
 817
 818   if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
 819   if(esdtrack) delete esdtrack;
 820   if(trackV0) delete trackV0;
 821
 822   return theCascade;
 823 }
 824
 825 //________________________________________________________________________
 826 AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::CallPrimaryVertex(AliAODv0 *v0, AliAODTrack *trk, AliAODEvent* aod)
 827 {
 828   //
 829   // Make an array of tracks which should not be used in primary vertex calculation and
 830   // Call PrimaryVertex function
 831   //
 832
 833   TObjArray *TrackArray = new TObjArray(3);
 834
 835   AliESDtrack *cptrk1 = new AliESDtrack((AliVTrack*)trk);
 836   TrackArray->AddAt(cptrk1,0);
 837
 838   AliESDtrack *cascptrack = new AliESDtrack((AliVTrack*)v0->GetDaughter(0));
 839   TrackArray->AddAt(cascptrack,1);
 840   AliESDtrack *cascntrack = new AliESDtrack((AliVTrack*)v0->GetDaughter(1));
 841   TrackArray->AddAt(cascntrack,2);
 842
 843   AliAODVertex *newvert  = PrimaryVertex(TrackArray,aod);
 844
 845   for(Int_t i=0;i<3;i++)
 846     {
 847       AliESDtrack *tesd = (AliESDtrack*)TrackArray->UncheckedAt(i);
 848       delete tesd;
 849     }
 850   TrackArray->Clear();
 851   delete TrackArray;
 852
 853   return newvert;
 854 }
 855
 856 //________________________________________________________________________
 857 AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::PrimaryVertex(const TObjArray *trkArray,
 858                                                                    AliVEvent *event)
 859 {
 860   //
 861   //Used only for pp
 862   //copied from AliAnalysisVertexingHF (except for the following 3 lines)
 863   //
 864
 865   Bool_t fRecoPrimVtxSkippingTrks = kTRUE;
 866   Bool_t fRmTrksFromPrimVtx = kFALSE;
 867
 868   AliESDVertex *vertexESD = 0;
 869   AliAODVertex *vertexAOD = 0;
 870
 871   //vertexESD = new AliESDVertex(*fV1);
 872
 873
 874   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
 875     // primary vertex from the input event
 876
 877     vertexESD = new AliESDVertex(*fV1);
 878
 879   } else {
 880     // primary vertex specific to this candidate
 881
 882     Int_t nTrks = trkArray->GetEntriesFast();
 883     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
 884
 885     if(fRecoPrimVtxSkippingTrks) {
 886       // recalculating the vertex
 887
 888       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
 889         Float_t diamondcovxy[3];
 890         event->GetDiamondCovXY(diamondcovxy);
 891         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
 892         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 893         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 894         vertexer->SetVtxStart(diamond);
 895         delete diamond; diamond=NULL;
 896         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
 897           vertexer->SetOnlyFitter();
 898       }
 899       Int_t skipped[1000];
 900       Int_t nTrksToSkip=0,id;
 901       AliExternalTrackParam *t = 0;
 902       for(Int_t i=0; i<nTrks; i++) {
 903         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
 904         id = (Int_t)t->GetID();
 905         if(id<0) continue;
 906         skipped[nTrksToSkip++] = id;
 907       }
 908       // TEMPORARY FIX
 909       // For AOD, skip also tracks without covariance matrix
 910       Double_t covtest[21];
 911       for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
 912         AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
 913         if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
 914           id = (Int_t)vtrack->GetID();
 915           if(id<0) continue;
 916           skipped[nTrksToSkip++] = id;
 917         }
 918       }
 919       for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
 920       //
 921       vertexer->SetSkipTracks(nTrksToSkip,skipped);
 922       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
 923
 924     } else if(fRmTrksFromPrimVtx && nTrks>0) {
 925       // removing the prongs tracks
 926
 927       TObjArray rmArray(nTrks);
 928       UShort_t *rmId = new UShort_t[nTrks];
 929       AliESDtrack *esdTrack = 0;
 930       AliESDtrack *t = 0;
 931       for(Int_t i=0; i<nTrks; i++) {
 932         t = (AliESDtrack*)trkArray->UncheckedAt(i);
 933         esdTrack = new AliESDtrack(*t);
 934         rmArray.AddLast(esdTrack);
 935         if(esdTrack->GetID()>=0) {
 936           rmId[i]=(UShort_t)esdTrack->GetID();
 937         } else {
 938           rmId[i]=9999;
 939         }
 940       }
 941       Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
 942       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
 943       delete [] rmId; rmId=NULL;
 944       rmArray.Delete();
 945
 946     }
 947
 948     delete vertexer; vertexer=NULL;
 949     if(!vertexESD) return vertexAOD;
 950     if(vertexESD->GetNContributors()<=0) {
 951       //AliDebug(2,"vertexing failed");
 952       delete vertexESD; vertexESD=NULL;
 953       return vertexAOD;
 954     }
 955
 956
 957   }
 958
 959   // convert to AliAODVertex
 960   Double_t pos[3],cov[6],chi2perNDF;
 961   vertexESD->GetXYZ(pos); // position
 962   vertexESD->GetCovMatrix(cov); //covariance matrix
 963   chi2perNDF = vertexESD->GetChi2toNDF();
 964   delete vertexESD; vertexESD=NULL;
 965
 966   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
 967
 968   return vertexAOD;
 969 }
 970
 971 //________________________________________________________________________
 972 AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::ReconstructSecondaryVertex(AliAODv0 *v0, AliAODTrack *part, AliAODEvent * aod)
 973 {
 974   //
 975   // Reconstruct secondary vertex from trkArray (Copied from AliAnalysisVertexingHF)
 976   //
 977   //------------------------------------------------
 978   // PrimaryVertex
 979   //------------------------------------------------
 980   AliAODVertex *primVertexAOD;
 981   Bool_t unsetvtx = kFALSE;
 982   if(fIspp){
 983     primVertexAOD = CallPrimaryVertex(v0,part,aod);
 984     if(!primVertexAOD){
 985       primVertexAOD = fVtx1;
 986     }else{
 987       unsetvtx = kTRUE;
 988     }
 989   }else{
 990     primVertexAOD = fVtx1;
 991   }
 992   if(!primVertexAOD) return 0x0;
 993
 994   //------------------------------------------------
 995   // Secondary vertex
 996   //------------------------------------------------
 997
 998   Double_t LcPx = part->Px()+v0->Px();
 999   Double_t LcPy = part->Py()+v0->Py();
1000   Double_t LcPt = TMath::Sqrt(LcPx*LcPx+LcPy*LcPy);
1001
1002   Double_t d0z0[2],covd0z0[3];
1003   part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1004   Double_t x0 = primVertexAOD->GetX();
1005   Double_t y0 = primVertexAOD->GetY();
1006   Double_t px0 = LcPx/LcPt;
1007   Double_t py0 = LcPy/LcPt;
1008   Double_t tx[3];
1009   part->GetXYZ(tx);
1010   Double_t x1 = tx[0];
1011   Double_t y1 = tx[1];
1012   part->GetPxPyPz(tx);
1013   Double_t px1 = tx[0];
1014   Double_t py1 = tx[1];
1015   Double_t pt1 = sqrt(px1*px1+py1*py1);
1016   px1 = px1/pt1;
1017   py1 = py1/pt1;
1018
1019   Double_t dx = x0 - x1;
1020   Double_t dy = y0 - y1;
1021
1022   Double_t Delta = -px0*py1+py0*px1;
1023   Double_t a0 = -9999.;
1024   if(Delta!=0)
1025     {
1026       a0 = 1./Delta * (py1 * dx - px1 * dy);
1027     }
1028   Double_t neovertx = x0 + a0 * px0;
1029   Double_t neoverty = y0 + a0 * py0;
1030   Double_t z0 = primVertexAOD->GetZ();
1031   Double_t neovertz = z0 + TMath::Abs(a0)*part->Pz()/part->Pt();
1032
1033   if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
1034
1035   Double_t pos[3],cov[6],chi2perNDF;
1036   pos[0]=neovertx;
1037   pos[1]=neoverty;
1038   pos[2]=neovertz;
1039   cov[0]=0.0;
1040   cov[1]=0.0;
1041   cov[2]=0.0;
1042   cov[3]=0.0;
1043   cov[4]=0.0;
1044   cov[5]=0.0;
1045   chi2perNDF=0.0;
1046   AliAODVertex *secVert = new AliAODVertex(pos,cov,chi2perNDF);
1047   if(!secVert){
1048     return 0x0;
1049   }
1050   return secVert;
1051 }