PWG4/PartCorr/AliV0Reader.cxx

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