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