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