PWG4/GammaConv/AliAnalysisTaskGammaConversion.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\r
  23 #include <TChain.h>\r
  24 \r
  25 // analysis\r
  26 #include "AliAnalysisTaskGammaConversion.h"\r
  27 #include "AliStack.h"\r
  28 #include "AliLog.h"\r
  29 #include <vector>\r
  30 \r
  31 class AliKFVertex;\r
  32 class AliAODHandler;\r
  33 class AliAODEvent;\r
  34 class ALiESDEvent;\r
  35 class AliMCEvent;\r
  36 class AliMCEventHandler;\r
  37 class AliESDInputHandler;\r
  38 class AliAnalysisManager;\r
  39 class Riostream;\r
  40 class TFile;\r
  41 class TInterpreter;\r
  42 class TSystem;\r
  43 class TROOT;\r
  44 \r
  45 ClassImp(AliAnalysisTaskGammaConversion)\r
  46 \r
  47 \r
  48 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():\r
  49   AliAnalysisTaskSE(),\r
  50   fV0Reader(NULL),\r
  51   fStack(NULL),\r
  52   fOutputContainer(NULL),\r
  53   fHistograms(NULL),\r
  54   fDoMCTruth(kFALSE),\r
  55   fMCAllGammas(),\r
  56   fMCPi0s(),\r
  57   fMCEtas(),\r
  58   fMCGammaChic(),\r
  59   fKFReconstructedGammas(),\r
  60   fElectronMass(-1),\r
  61   fGammaMass(-1),\r
  62   fPi0Mass(-1),\r
  63   fEtaMass(-1),\r
  64   fGammaWidth(-1),\r
  65   fPi0Width(-1),\r
  66   fEtaWidth(-1),\r
  67   fCalculateBackground(kFALSE)\r
  68 {\r
  69   // Default constructor\r
  70   // Common I/O in slot 0\r
  71   DefineInput (0, TChain::Class());\r
  72   DefineOutput(0, TTree::Class());\r
  73 \r
  74   // Your private output\r
  75   DefineOutput(1, TList::Class());\r
  76 }\r
  77 \r
  78 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
  79   AliAnalysisTaskSE(name),\r
  80   fV0Reader(NULL),\r
  81   fStack(NULL),\r
  82   fOutputContainer(0x0),\r
  83   fHistograms(NULL),\r
  84   fDoMCTruth(kFALSE),\r
  85   fMCAllGammas(),\r
  86   fMCPi0s(),\r
  87   fMCEtas(),\r
  88   fMCGammaChic(),\r
  89   fKFReconstructedGammas(),\r
  90   fElectronMass(-1),\r
  91   fGammaMass(-1),\r
  92   fPi0Mass(-1),\r
  93   fEtaMass(-1),\r
  94   fGammaWidth(-1),\r
  95   fPi0Width(-1),\r
  96   fEtaWidth(-1),\r
  97   fCalculateBackground(kFALSE)\r
  98 {\r
  99   // Common I/O in slot 0\r
 100   DefineInput (0, TChain::Class());\r
 101   DefineOutput(0, TTree::Class());\r
 102   \r
 103   // Your private output\r
 104   DefineOutput(1, TList::Class());\r
 105 }\r
 106 \r
 107 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
 108 {\r
 109   // Remove all pointers\r
 110  \r
 111   if(fOutputContainer){\r
 112     fOutputContainer->Clear() ; \r
 113     delete fOutputContainer ;\r
 114   }\r
 115   if(fHistograms){\r
 116     delete fHistograms;\r
 117   }\r
 118   if(fV0Reader){\r
 119     delete fV0Reader;\r
 120   }\r
 121 }\r
 122 \r
 123 \r
 124 void AliAnalysisTaskGammaConversion::Init()\r
 125 {\r
 126   // Initialization\r
 127   AliLog::SetGlobalLogLevel(AliLog::kError);\r
 128 }\r
 129 \r
 130 \r
 131 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
 132 {\r
 133   // Execute analysis for current event\r
 134   \r
 135   ConnectInputData("");\r
 136   \r
 137   //clear vectors\r
 138   fMCAllGammas.clear();\r
 139   fMCPi0s.clear();\r
 140   fMCEtas.clear();\r
 141   fMCGammaChic.clear();\r
 142 \r
 143   fKFReconstructedGammas.clear();\r
 144 \r
 145   //Clear the data in the v0Reader\r
 146   fV0Reader->UpdateEventByEventData();\r
 147 \r
 148   // Process the MC information\r
 149   if(fDoMCTruth){\r
 150     ProcessMCData();\r
 151   }\r
 152 \r
 153   // Process the v0 information\r
 154   ProcessV0s();\r
 155 \r
 156   //calculate background if flag is set\r
 157   if(fCalculateBackground){\r
 158     CalculateBackground();\r
 159   }\r
 160 \r
 161   // Process reconstructed gammas\r
 162   ProcessGammasForNeutralMesonAnalysis();\r
 163 \r
 164   PostData(1, fOutputContainer);\r
 165   \r
 166 }\r
 167 \r
 168 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
 169   // see header file for documentation\r
 170 \r
 171   if(fV0Reader == NULL){\r
 172     // Write warning here cuts and so on are default if this ever happens\r
 173   }\r
 174   fV0Reader->Initialize();\r
 175 }\r
 176 \r
 177 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
 178   // see header file for documentation\r
 179   \r
 180   fStack = fV0Reader->GetMCStack();\r
 181 \r
 182   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
 183     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
 184 \r
 185     if (!particle) {\r
 186       //print warning here\r
 187       continue;\r
 188     }\r
 189     \r
 190     if(particle->Pt()<fV0Reader->GetPtCut()){\r
 191       continue;\r
 192     }\r
 193     \r
 194     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut()){\r
 195       continue;\r
 196     }\r
 197 \r
 198     if(particle->R()>fV0Reader->GetMaxRCut()){ // cuts on distance from collision point\r
 199       continue;\r
 200     }\r
 201 \r
 202     Double_t tmpPhi=particle->Phi();\r
 203     if(particle->Phi()> TMath::Pi()){\r
 204       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
 205     }\r
 206 \r
 207     \r
 208     //process the gammas\r
 209     if (particle->GetPdgCode()== 22){\r
 210       fMCAllGammas.push_back(particle);\r
 211       if(particle->GetMother(0)>-1){ //Means we have a mother\r
 212         if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 ){//Checks for a non gamma mother.\r
 213           fHistograms->FillHistogram("MC_Gamma_Energy", particle->Energy());\r
 214           fHistograms->FillHistogram("MC_Gamma_Pt", particle->Pt());\r
 215           fHistograms->FillHistogram("MC_Gamma_Eta", particle->Eta());\r
 216           \r
 217           fHistograms->FillHistogram("MC_Gamma_Phi", tmpPhi);\r
 218 \r
 219           //adding the conversion points from all gammas with e+e- daughters\r
 220           if(particle->GetNDaughters() >= 2){\r
 221             TParticle* daughter0 = NULL;\r
 222             TParticle* daughter1 = NULL;\r
 223             \r
 224             for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
 225               TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
 226               if(tmpDaughter->GetUniqueID() == 5){\r
 227                 if(tmpDaughter->GetPdgCode() == 11){\r
 228                   daughter0 = tmpDaughter;\r
 229                 }\r
 230                 else if(tmpDaughter->GetPdgCode() == -11){\r
 231                   daughter1 = tmpDaughter;\r
 232                 }\r
 233               }\r
 234             }\r
 235 \r
 236             if(daughter0 == NULL || daughter1 == NULL){ // means we do not have two daughters from pair production\r
 237               continue;\r
 238             }\r
 239 \r
 240             if(daughter0->R()>fV0Reader->GetMaxRCut() || daughter1->R()>fV0Reader->GetMaxRCut()){\r
 241               continue;\r
 242             }\r
 243             \r
 244             if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||\r
 245                (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
 246 \r
 247               // begin Mapping \r
 248               Int_t rBin    = fHistograms->GetRBin(daughter0->R());\r
 249               Int_t phiBin  = fHistograms->GetPhiBin(daughter0->Phi());\r
 250               \r
 251               TString nameMCMappingPhiR="";\r
 252               nameMCMappingPhiR.Form("MC_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 253               fHistograms->FillHistogram(nameMCMappingPhiR, daughter0->Vz(), particle->Eta());\r
 254               \r
 255               TString nameMCMappingPhi="";\r
 256               nameMCMappingPhi.Form("MC_EP_Mapping-Phi%02d",phiBin);\r
 257               fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
 258               \r
 259               TString nameMCMappingR="";\r
 260               nameMCMappingR.Form("MC_EP_Mapping-R%02d",rBin);\r
 261               fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
 262               \r
 263               TString nameMCMappingPhiInR="";\r
 264               nameMCMappingPhiInR.Form("MC_EP_Mapping_Phi_R-%02d",rBin);\r
 265               fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
 266               //end mapping\r
 267 \r
 268               fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
 269               fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
 270               fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
 271               fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
 272             }// end if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||....... approx 20 lines above\r
 273           }// end if(particle->GetNDaughters() >= 2){\r
 274         } // end if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 )\r
 275         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
 276             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
 277             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
 278             ){ \r
 279           fMCGammaChic.push_back(particle);\r
 280         }\r
 281       }// end if(particle->GetMother(0)>-1)\r
 282       else{//means we have a primary particle\r
 283         fHistograms->FillHistogram("MC_DirectGamma_Energy",particle->Energy());\r
 284         fHistograms->FillHistogram("MC_DirectGamma_Pt", particle->Pt());\r
 285         fHistograms->FillHistogram("MC_DirectGamma_Eta", particle->Eta());\r
 286         fHistograms->FillHistogram("MCDirectGammaPhi", tmpPhi);\r
 287 \r
 288         //adding the conversion points from all gammas with e+e- daughters\r
 289         if(particle->GetNDaughters() == 2){\r
 290           TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
 291           TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
 292           if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11) ||\r
 293              (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
 294             \r
 295             fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
 296             fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
 297             fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
 298             fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
 299 \r
 300           }\r
 301         }\r
 302       }// end else\r
 303     }// end if (particle->GetPdgCode()== 22){\r
 304     else if (TMath::Abs(particle->GetPdgCode())== 11){ // Means we have an electron or a positron\r
 305       if(particle->GetMother(0)>-1){ // means we have a mother\r
 306         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==22 ){ // Means we have a gamma mother\r
 307           if(particle->GetPdgCode() == 11){//electron \r
 308             fHistograms->FillHistogram("MC_E_Energy", particle->Energy());\r
 309             fHistograms->FillHistogram("MC_E_Pt", particle->Pt());\r
 310             fHistograms->FillHistogram("MC_E_Eta", particle->Eta());\r
 311             fHistograms->FillHistogram("MC_E_Phi", tmpPhi);\r
 312           }\r
 313           if(particle->GetPdgCode() == -11){//positron \r
 314             fHistograms->FillHistogram("MC_P_Energy", particle->Energy());\r
 315             fHistograms->FillHistogram("MC_P_Pt", particle->Pt());\r
 316             fHistograms->FillHistogram("MC_P_Eta", particle->Eta());\r
 317             fHistograms->FillHistogram("MC_P_Phi", tmpPhi);\r
 318           }\r
 319         }\r
 320       }\r
 321     } // end else if (TMath::Abs(particle->GetPdgCode())== 11)\r
 322     else if(particle->GetNDaughters() == 2){\r
 323 \r
 324       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
 325       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
 326       if(daughter0->GetPdgCode() == 22 && daughter1->GetPdgCode() == 22){//check for gamma gamma daughters\r
 327         \r
 328         if(particle->GetPdgCode()==111){//Pi0\r
 329           \r
 330           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
 331 \r
 332           if( iTracks >= fStack->GetNprimary()){\r
 333             \r
 334             fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
 335 \r
 336             fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
 337             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
 338             fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
 339             fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
 340             fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
 341             fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 342             fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
 343           }\r
 344           else{\r
 345             fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());\r
 346 \r
 347             fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
 348             fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
 349             fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
 350             fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
 351             fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
 352             fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 353             fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 354           }\r
 355         }\r
 356         else if(particle->GetPdgCode()==221){//Eta\r
 357           fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
 358 \r
 359           fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
 360           fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
 361           fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
 362           fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
 363           fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
 364           fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 365           fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 366         }\r
 367         \r
 368         //the match data should be filled no matter which mother the gamma-gamma comes from\r
 369         fHistograms->FillHistogram("MC_Match_Gamma_R", particle->R());\r
 370         fHistograms->FillHistogram("MC_Match_Gamma_ZR", particle->Vz(),particle->R());\r
 371         fHistograms->FillHistogram("MC_Match_Gamma_XY", particle->Vx(),particle->Vy());\r
 372         fHistograms->FillHistogram("MC_Match_Gamma_Mass", particle->GetCalcMass());\r
 373         fHistograms->FillHistogram("MC_Match_Gamma_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 374         fHistograms->FillHistogram("MC_Match_Gamma_Energy", particle->Energy());\r
 375         fHistograms->FillHistogram("MC_Match_Gamma_Pt", particle->Pt());\r
 376         fHistograms->FillHistogram("MC_Match_Gamma_Eta", particle->Eta());\r
 377         fHistograms->FillHistogram("MC_Match_Gamma_Phi",tmpPhi);\r
 378       }\r
 379     }// end else if(particle->GetNDaughters() == 2)\r
 380   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
 381 } // end ProcessMCData\r
 382 \r
 383 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
 384   // see header file for documentation\r
 385 \r
 386   Int_t nSurvivingV0s=0;\r
 387   while(fV0Reader->NextV0()){\r
 388     nSurvivingV0s++;\r
 389     //-------------------------- filling v0 information -------------------------------------\r
 390     fHistograms->FillHistogram("ESD_EP_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
 391     fHistograms->FillHistogram("ESD_EP_R", fV0Reader->GetXYRadius());\r
 392     fHistograms->FillHistogram("ESD_EP_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 393     fHistograms->FillHistogram("ESD_EP_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 394     \r
 395     \r
 396     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
 397     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
 398     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
 399     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
 400     \r
 401     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
 402     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
 403     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
 404     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
 405     \r
 406     fHistograms->FillHistogram("ESD_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 407     fHistograms->FillHistogram("ESD_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 408     fHistograms->FillHistogram("ESD_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
 409     fHistograms->FillHistogram("ESD_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 410 \r
 411 \r
 412     // begin mapping\r
 413     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
 414     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
 415     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
 416     \r
 417     TString nameESDMappingPhiR="";\r
 418     nameESDMappingPhiR.Form("ESD_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 419     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
 420 \r
 421     TString nameESDMappingPhi="";\r
 422     nameESDMappingPhi.Form("ESD_EP_Mapping-Phi%02d",phiBin);\r
 423     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
 424 \r
 425     TString nameESDMappingR="";\r
 426     nameESDMappingR.Form("ESD_EP_Mapping-R%02d",rBin);\r
 427     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
 428 \r
 429     TString nameESDMappingPhiInR="";\r
 430     nameESDMappingPhiInR.Form("ESD_EP_Mapping_Phi_R-%02d",rBin);\r
 431     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
 432     // end mapping\r
 433     \r
 434     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
 435 \r
 436     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
 437     if(fDoMCTruth){\r
 438       if(fV0Reader->HasSameMCMother() == kFALSE){\r
 439         continue;\r
 440       }\r
 441  \r
 442       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
 443       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
 444 \r
 445       if(negativeMC->GetPdgCode()!=11 || positiveMC->GetPdgCode()!=-11){\r
 446         continue;\r
 447       }\r
 448       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
 449         fHistograms->FillHistogram("ESD_Match_Gamma_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 450         fHistograms->FillHistogram("ESD_Match_Gamma_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 451         fHistograms->FillHistogram("ESD_Match_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 452         fHistograms->FillHistogram("ESD_Match_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 453         fHistograms->FillHistogram("ESD_Match_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
 454 \r
 455         fHistograms->FillHistogram("ESD_Match_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 456         fHistograms->FillHistogram("ESD_Match_Gamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 457         fHistograms->FillHistogram("ESD_Match_Gamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 458         fHistograms->FillHistogram("ESD_Match_Gamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 459         fHistograms->FillHistogram("ESD_Match_Gamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 460         fHistograms->FillHistogram("ESD_Match_Gamma_R", fV0Reader->GetXYRadius());\r
 461         fHistograms->FillHistogram("ESD_Match_Gamma_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 462 \r
 463         //resolution\r
 464         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
 465         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
 466         Double_t resdPt = 0;\r
 467         if(mcpt != 0){\r
 468           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
 469         }\r
 470 \r
 471         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
 472         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
 473         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
 474         \r
 475         Double_t resdZ = 0;\r
 476         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
 477           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
 478         }\r
 479         \r
 480         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
 481         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
 482         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
 483         \r
 484         Double_t resdR = 0;\r
 485         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
 486           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
 487         }\r
 488         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
 489         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
 490         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
 491         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
 492       }\r
 493     }\r
 494   }\r
 495   fHistograms->FillHistogram("NumberOfSurvivingV0s", nSurvivingV0s);\r
 496   fHistograms->FillHistogram("NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
 497 }\r
 498 \r
 499 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
 500   // see header file for documentation\r
 501 \r
 502   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
 503     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
 504       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
 505       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
 506 \r
 507       \r
 508       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
 509 \r
 510       Double_t massTwoGammaCandidate =0.;\r
 511       Double_t widthTwoGammaCandidate = 0.;\r
 512       Double_t chi2TwoGammaCandidate =10000.;   \r
 513       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
 514       if(twoGammaCandidate->GetNDF()>0){\r
 515         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
 516         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){\r
 517 \r
 518           TVector3 vectorTwoGammaCandidate(twoGammaCandidate->Px(),twoGammaCandidate->Py(),twoGammaCandidate->Pz());\r
 519 \r
 520           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);\r
 521 \r
 522           //Calculating by hand the radius\r
 523           Double_t tmpX= twoGammaCandidate->GetX();\r
 524           Double_t tmpY= twoGammaCandidate->GetY();\r
 525           \r
 526           Double_t radiusTwoGammaCandidate = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
 527 \r
 528           fHistograms->FillHistogram("ESD_TwoGammaCombination_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
 529           fHistograms->FillHistogram("ESD_TwoGammaCombination_Energy", twoGammaCandidate->GetE());\r
 530           fHistograms->FillHistogram("ESD_TwoGammaCombination_Pt", sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
 531           fHistograms->FillHistogram("ESD_TwoGammaCombination_Eta", vectorTwoGammaCandidate.Eta());\r
 532           fHistograms->FillHistogram("ESD_TwoGammaCombination_Phi", vectorTwoGammaCandidate.Phi());\r
 533           fHistograms->FillHistogram("ESD_TwoGammaCombination_Mass", massTwoGammaCandidate);\r
 534           fHistograms->FillHistogram("ESD_TwoGammaCombination_R", radiusTwoGammaCandidate);\r
 535           fHistograms->FillHistogram("ESD_TwoGammaCombination_ZR", tmpY, radiusTwoGammaCandidate);\r
 536           fHistograms->FillHistogram("ESD_TwoGammaCombination_XY", tmpX, tmpY);\r
 537           fHistograms->FillHistogram("InvMass_vs_Pt_Spectra",massTwoGammaCandidate ,sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
 538         }\r
 539       }\r
 540       delete twoGammaCandidate;\r
 541     }\r
 542   }\r
 543 }\r
 544 \r
 545 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
 546   // see header file for documentation\r
 547 \r
 548   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
 549   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
 550   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
 551     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
 552     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
 553       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
 554 \r
 555       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
 556 \r
 557       Double_t massBG =0.;\r
 558       Double_t widthBG = 0.;\r
 559       Double_t chi2BG =10000.;  \r
 560       backgroundCandidate->GetMass(massBG,widthBG);\r
 561       if(backgroundCandidate->GetNDF()>0){\r
 562         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
 563         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
 564 \r
 565           TVector3 vectorBGCandidate(backgroundCandidate->Px(),backgroundCandidate->Py(),backgroundCandidate->Pz());\r
 566 \r
 567           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
 568 \r
 569           //Calculating by hand the radius (find a better way)\r
 570           Double_t tmpX= backgroundCandidate->GetX();\r
 571           Double_t tmpY= backgroundCandidate->GetY();\r
 572           \r
 573           Double_t radiusBG = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
 574 \r
 575           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
 576           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
 577           fHistograms->FillHistogram("ESD_Background_Pt", sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
 578           fHistograms->FillHistogram("ESD_Background_Eta", vectorBGCandidate.Eta());\r
 579           fHistograms->FillHistogram("ESD_Background_Phi", vectorBGCandidate.Phi());\r
 580           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
 581           fHistograms->FillHistogram("ESD_Background_R", radiusBG);\r
 582           fHistograms->FillHistogram("ESD_Background_ZR", tmpY, radiusBG);\r
 583           fHistograms->FillHistogram("ESD_Background_XY", tmpX, tmpY);\r
 584           fHistograms->FillHistogram("Background_InvMass_vs_Pt_Spectra",massBG,sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
 585         }\r
 586       }\r
 587       delete backgroundCandidate;   \r
 588     }\r
 589   }\r
 590 }\r
 591 \r
 592 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
 593 {\r
 594   // Terminate analysis\r
 595   //\r
 596   AliDebug(1,"Do nothing in Terminate");\r
 597 }\r
 598 \r
 599 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
 600 {\r
 601   // Create the output container\r
 602   if(fOutputContainer != NULL){\r
 603     delete fOutputContainer;\r
 604     fOutputContainer = NULL;\r
 605   }\r
 606   if(fOutputContainer == NULL){\r
 607     fOutputContainer = new TList();\r
 608   }\r
 609   fHistograms->GetOutputContainer(fOutputContainer);\r
 610   fOutputContainer->SetName(GetName());\r
 611 }\r
 612 \r
 613 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const{\r
 614   //helper function\r
 615   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
 616   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
 617   return v3D0.Angle(v3D1);\r
 618 }\r