PWG4/AliV0Reader.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *
   5  * Version 1.0                                                            *
   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 /**
  16  * Class for reading V0's
  17  */
  18
  19 // --- ROOT system ---
  20 #include <TFormula.h>
  21 #include <TMath.h>
  22
  23 //---- ANALYSIS system ----
  24 #include "AliV0Reader.h"
  25 #include "AliAnalysisManager.h"
  26 #include "AliESDInputHandler.h"
  27 #include "AliESDv0.h"
  28 #include "AliMCEvent.h"
  29 #include "AliKFVertex.h"
  30 #include <iostream>
  31
  32 using namespace std;
  33
  34 ClassImp(AliV0Reader)
  35
  36
  37
  38 AliV0Reader::AliV0Reader() :
  39   TObject(),
  40   fCurrentEventGoodV0s(),
  41   fPreviousEventGoodV0s(),
  42   fMCStack(NULL),
  43   fMCTruth(NULL),
  44   fChain(NULL),
  45   fESDHandler(NULL),
  46   fESDEvent(NULL),
  47   fHistograms(NULL),
  48   fCurrentV0IndexNumber(0),
  49   fCurrentV0(NULL),
  50   fCurrentNegativeKFParticle(NULL),
  51   fCurrentPositiveKFParticle(NULL),
  52   fCurrentMotherKFCandidate(NULL),
  53   fCurrentNegativeESDTrack(NULL),
  54   fCurrentPositiveESDTrack(NULL),
  55   fNegativeTrackLorentzVector(NULL),
  56   fPositiveTrackLorentzVector(NULL),
  57   fMotherCandidateLorentzVector(NULL),
  58   fCurrentXValue(0),
  59   fCurrentYValue(0),
  60   fCurrentZValue(0),
  61   fPositiveTrackPID(0),
  62   fNegativeTrackPID(0),
  63   fNegativeMCParticle(NULL),
  64   fPositiveMCParticle(NULL),
  65   fMotherMCParticle(NULL),
  66   fMotherCandidateKFMass(0),
  67   fMotherCandidateKFWidth(0),
  68   fUseKFParticle(kTRUE),
  69   fUseESDTrack(kFALSE),
  70   fDoMC(kFALSE),
  71   fMaxR(10000),// 100 meter(outside of ALICE)
  72   fEtaCut(0.),
  73   fPtCut(0.),
  74   fChi2Cut(0.),
  75   fPIDProbabilityCutNegativeParticle(0),
  76   fPIDProbabilityCutPositiveParticle(0),
  77   fXVertexCut(0.),
  78   fYVertexCut(0.),
  79   fZVertexCut(0.),
  80   fNSigmaMass(0.),
  81   fUseImprovedVertex(kFALSE)
  82 {
  83
  84 }
  85
  86
  87 AliV0Reader::AliV0Reader(const AliV0Reader & original) :
  88   TObject(original),
  89   fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),
  90   fPreviousEventGoodV0s(original.fPreviousEventGoodV0s),
  91   fMCStack(original.fMCStack),
  92   fMCTruth(original.fMCTruth),
  93   fChain(original.fChain),
  94   fESDHandler(original.fESDHandler),
  95   fESDEvent(original.fESDEvent),
  96   fHistograms(original.fHistograms),
  97   fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),
  98   fCurrentV0(original.fCurrentV0),
  99   fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),
 100   fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),
 101   fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),
 102   fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),
 103   fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),
 104   fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),
 105   fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),
 106   fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),
 107   fCurrentXValue(original.fCurrentXValue),
 108   fCurrentYValue(original.fCurrentYValue),
 109   fCurrentZValue(original.fCurrentZValue),
 110   fPositiveTrackPID(original.fPositiveTrackPID),
 111   fNegativeTrackPID(original.fNegativeTrackPID),
 112   fNegativeMCParticle(original.fNegativeMCParticle),
 113   fPositiveMCParticle(original.fPositiveMCParticle),
 114   fMotherMCParticle(original.fMotherMCParticle),
 115   fMotherCandidateKFMass(original.fMotherCandidateKFMass),
 116   fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),
 117   fUseKFParticle(kTRUE),
 118   fUseESDTrack(kFALSE),
 119   fDoMC(kFALSE),
 120   fMaxR(original.fMaxR),
 121   fEtaCut(original.fEtaCut),
 122   fPtCut(original.fPtCut),
 123   fChi2Cut(original.fChi2Cut),
 124   fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),
 125   fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),
 126   fXVertexCut(original.fXVertexCut),
 127   fYVertexCut(original.fYVertexCut),
 128   fZVertexCut(original.fZVertexCut),
 129   fNSigmaMass(original.fNSigmaMass),
 130   fUseImprovedVertex(original.fUseImprovedVertex)
 131 {
 132
 133 }
 134
 135
 136 AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)
 137 {
 138   // assignment operator
 139   return *this;
 140 }
 141
 142 void AliV0Reader::Initialize(){
 143   // Get the input handler from the manager
 144   fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 145   if(fESDHandler == NULL){
 146     //print warning here
 147   }
 148
 149   // Get pointer to esd event from input handler
 150   fESDEvent = fESDHandler->GetEvent();
 151   if(fESDEvent == NULL){
 152     //print warning here
 153   }
 154
 155   //Get pointer to MCTruth
 156   fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 157   if(fMCTruth == NULL){
 158     //print warning here
 159   }
 160
 161   //Get pointer to the mc stack
 162   fMCStack = fMCTruth->MCEvent()->Stack();
 163   if(fMCStack == NULL){
 164     //print warning here
 165   }
 166
 167   AliKFParticle::SetField(fESDEvent->GetMagneticField());
 168
 169 }
 170
 171 AliESDv0* AliV0Reader::GetV0(Int_t index){
 172   fCurrentV0 = fESDEvent->GetV0(index);
 173   UpdateV0Information();
 174   return fCurrentV0;
 175 }
 176
 177
 178 Bool_t AliV0Reader::NextV0(){
 179   Bool_t iResult=kFALSE;
 180   while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){
 181     fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);
 182
 183    //checks if on the fly mode is set
 184     if ( !fCurrentV0->GetOnFlyStatus() ){
 185       fCurrentV0IndexNumber++;
 186       if(fHistograms->fV0MassDebugCut1){fHistograms->fV0MassDebugCut1->Fill(GetMotherCandidateMass());}
 187       continue;
 188     }
 189
 190     if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {//checks if we have a vertex
 191       fCurrentV0IndexNumber++;
 192       if(fHistograms->fV0MassDebugCut2){fHistograms->fV0MassDebugCut2->Fill(GetMotherCandidateMass());}
 193       continue;
 194     }
 195
 196     if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){
 197       fCurrentV0IndexNumber++;
 198       if(fHistograms->fV0MassDebugCut3){fHistograms->fV0MassDebugCut3->Fill(GetMotherCandidateMass());}
 199       continue;
 200     }
 201
 202
 203     fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
 204
 205     if(GetXYRadius()>fMaxR){ // cuts on distance from collision point
 206       fCurrentV0IndexNumber++;
 207       if(fHistograms->fV0MassDebugCut4){fHistograms->fV0MassDebugCut4->Fill(GetMotherCandidateMass());}
 208       continue;
 209     }
 210
 211     UpdateV0Information();
 212
 213     if(fUseKFParticle){
 214       if(fCurrentMotherKFCandidate->GetNDF()<=0){
 215         fCurrentV0IndexNumber++;
 216       if(fHistograms->fV0MassDebugCut5){fHistograms->fV0MassDebugCut5->Fill(GetMotherCandidateMass());}
 217         continue;
 218       }
 219       Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();
 220       if(chi2V0 > fChi2Cut || chi2V0 <=0){
 221         fCurrentV0IndexNumber++;
 222       if(fHistograms->fV0MassDebugCut6){fHistograms->fV0MassDebugCut6->Fill(GetMotherCandidateMass());}
 223         continue;
 224       }
 225
 226       if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){
 227         fCurrentV0IndexNumber++;
 228       if(fHistograms->fV0MassDebugCut7){fHistograms->fV0MassDebugCut7->Fill(GetMotherCandidateMass());}
 229         continue;
 230       }
 231
 232       if(fMotherCandidateLorentzVector->Pt()<fPtCut){
 233         fCurrentV0IndexNumber++;
 234       if(fHistograms->fV0MassDebugCut8){fHistograms->fV0MassDebugCut8->Fill(GetMotherCandidateMass());}
 235         continue;
 236       }
 237
 238     }
 239     else if(fUseESDTrack){
 240       //TODO
 241     }
 242
 243     iResult=kTRUE;//means we have a v0 who survived all the cuts applied
 244
 245     fCurrentV0IndexNumber++;
 246
 247     break;
 248   }
 249   return iResult;
 250 }
 251
 252 void AliV0Reader::UpdateV0Information(){
 253     if(fCurrentNegativeKFParticle != NULL){
 254       delete fCurrentNegativeKFParticle;
 255       //      fCurrentNegativeKFParticle = NULL;
 256     }
 257     fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);
 258
 259     if(fCurrentPositiveKFParticle != NULL){
 260       delete fCurrentPositiveKFParticle;
 261       //      fCurrentPositiveKFParticle = NULL;
 262     }
 263     fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);
 264
 265     if(fCurrentMotherKFCandidate != NULL){
 266       //      cout<<"fCurrentMotherKFCandidate: "<<fCurrentMotherKFCandidate<<endl;
 267       delete fCurrentMotherKFCandidate;
 268       //      fCurrentMotherKFCandidate = NULL;
 269     }
 270     fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
 271
 272     fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
 273
 274     fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
 275
 276     if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
 277       fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
 278     }
 279
 280     if(fUseImprovedVertex == kTRUE){
 281       AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
 282       primaryVertexImproved+=*fCurrentMotherKFCandidate;
 283       fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
 284     }
 285
 286     fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);
 287
 288
 289     if(fNegativeTrackLorentzVector != NULL){
 290       delete fNegativeTrackLorentzVector;
 291     }
 292     if(fUseKFParticle){
 293       fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());
 294     }
 295     else if(fUseESDTrack){
 296       fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());
 297     }
 298
 299     if(fPositiveTrackLorentzVector != NULL){
 300       delete fPositiveTrackLorentzVector;
 301     }
 302     if(fUseKFParticle){
 303       fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());
 304     }
 305     else if(fUseESDTrack){
 306       fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());
 307     }
 308
 309     if(fMotherCandidateLorentzVector != NULL){
 310       delete fMotherCandidateLorentzVector;
 311     }
 312     if(fUseKFParticle){
 313       fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
 314         // new TLorentzVector(fCurrentMotherKFCandidate->Px(),fCurrentMotherKFCandidate->Py(),fCurrentMotherKFCandidate->Pz());
 315     }
 316     else if(fUseESDTrack){
 317       fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
 318         // new TLorentzVector(fCurrentNegativeESDTrack->Px()+fCurrentPositiveESDTrack->Px(),fCurrentNegativeESDTrack->Py()+fCurrentPositiveESDTrack->Py(),fCurrentNegativeESDTrack->Pz()+fCurrentPositiveESDTrack->Pz());
 319     }
 320
 321     if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
 322       fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.);
 323     }
 324
 325     if(fDoMC){
 326       fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
 327       fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
 328     }
 329 }
 330
 331 Bool_t AliV0Reader::HasSameMCMother(){
 332   Bool_t iResult = kFALSE;
 333   if(fDoMC == kTRUE){
 334     if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){
 335       if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))
 336         fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));
 337         iResult = kTRUE;
 338       }
 339   }
 340   return iResult;
 341 }
 342
 343 AliKFParticle* AliV0Reader::GetNegativeKFParticle(){
 344   return fCurrentNegativeKFParticle;
 345 }
 346
 347 AliKFParticle* AliV0Reader::GetPositiveKFParticle(){
 348   return fCurrentPositiveKFParticle;
 349 }
 350
 351 AliKFParticle* AliV0Reader::GetMotherCandidateKFCombination(){
 352   return fCurrentMotherKFCandidate;
 353 }
 354
 355 Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){
 356   Bool_t iResult=kFALSE;
 357
 358   Double_t *posProbArray = new Double_t[10];
 359   Double_t *negProbArray = new Double_t[10];
 360   AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());
 361   AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());
 362
 363   negTrack->GetTPCpid(negProbArray);
 364   posTrack->GetTPCpid(posProbArray);
 365
 366   if(negProbArray!=NULL && posProbArray!=NULL){
 367     if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){
 368       iResult=kTRUE;
 369     }
 370   }
 371   delete [] posProbArray;
 372   delete [] negProbArray;
 373   return iResult;
 374 }
 375
 376 void AliV0Reader::UpdateEventByEventData(){
 377   fPreviousEventGoodV0s.clear();
 378   fPreviousEventGoodV0s = fCurrentEventGoodV0s;
 379   fCurrentEventGoodV0s.clear();
 380
 381   fCurrentV0IndexNumber=0;
 382 }
 383
 384 Double_t AliV0Reader::GetNegativeTrackPhi(){
 385   Double_t offset=0;
 386   if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
 387     offset = -2*TMath::Pi();
 388   }
 389   return fNegativeTrackLorentzVector->Phi()+offset;
 390 }
 391
 392 Double_t AliV0Reader::GetPositiveTrackPhi(){
 393   Double_t offset=0;
 394   if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
 395     offset = -2*TMath::Pi();
 396   }
 397   return fPositiveTrackLorentzVector->Phi()+offset;
 398 }
 399
 400 Double_t AliV0Reader::GetMotherCandidatePhi(){
 401   Double_t offset=0;
 402   if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
 403     offset = -2*TMath::Pi();
 404   }
 405   return fMotherCandidateLorentzVector->Phi()+offset;
 406 }
 407
 408 Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
 409
 410   Int_t iResult = 10; // Unknown particle
 411
 412   if(chargeOfTrack==-1){ //negative track
 413     switch(abs(fNegativeTrackPID)){
 414     case 11:       //electron
 415       iResult = 0;
 416       break;
 417     case 13:       //muon
 418       iResult = 1;
 419       break;
 420     case 211:      //pion
 421       iResult = 2;
 422       break;
 423     case 321:      //kaon
 424       iResult = 3;
 425       break;
 426     case 2212:     //proton
 427       iResult = 4;
 428       break;
 429     case 22:       //photon
 430       iResult = 5;
 431       break;
 432     case 111:      //pi0
 433       iResult = 6;
 434       break;
 435     case 2112:     //neutron
 436       iResult = 7;
 437       break;
 438     case 311:      //K0
 439       iResult = 8;
 440       break;
 441
 442       //Put in here for kSPECIES::kEleCon  ????
 443     }
 444   }
 445   else if(chargeOfTrack==1){ //positive track
 446     switch(abs(fPositiveTrackPID)){
 447     case 11:       //electron
 448       iResult = 0;
 449       break;
 450     case 13:       //muon
 451       iResult = 1;
 452       break;
 453     case 211:      //pion
 454       iResult = 2;
 455       break;
 456     case 321:      //kaon
 457       iResult = 3;
 458       break;
 459     case 2212:     //proton
 460       iResult = 4;
 461       break;
 462     case 22:       //photon
 463       iResult = 5;
 464       break;
 465     case 111:      //pi0
 466       iResult = 6;
 467       break;
 468     case 2112:     //neutron
 469       iResult = 7;
 470       break;
 471     case 311:      //K0
 472       iResult = 8;
 473       break;
 474
 475       //Put in here for kSPECIES::kEleCon  ????
 476     }
 477   }
 478   else{
 479     //Wrong parameter.. Print warning
 480   }
 481   return iResult;
 482 }