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 "TNtuple.h"\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   fESDEvent(NULL),      \r
  53   fOutputContainer(NULL),\r
  54   fHistograms(NULL),\r
  55   fDoMCTruth(kFALSE),\r
  56   fMCAllGammas(),\r
  57   fMCPi0s(),\r
  58   fMCEtas(),\r
  59   fMCGammaChic(),\r
  60   fKFReconstructedGammas(),\r
  61   fIsTrueReconstructedGammas(),\r
  62   fElectronv1(),\r
  63   fElectronv2(),\r
  64   fCurrentEventPosElectron(),\r
  65   fPreviousEventPosElectron(),\r
  66   fCurrentEventNegElectron(),\r
  67   fPreviousEventNegElectron(),\r
  68   fKFReconstructedGammasCut(),                  \r
  69   fPreviousEventTLVNegElectron(),\r
  70   fPreviousEventTLVPosElectron(),       \r
  71   fElectronMass(-1),\r
  72   fGammaMass(-1),\r
  73   fPi0Mass(-1),\r
  74   fEtaMass(-1),\r
  75   fGammaWidth(-1),\r
  76   fPi0Width(-1),\r
  77   fEtaWidth(-1),\r
  78   fMinOpeningAngleGhostCut(0.),\r
  79   fCalculateBackground(kFALSE),\r
  80   fWriteNtuple(kFALSE),\r
  81   fGammaNtuple(NULL),\r
  82   fNeutralMesonNtuple(NULL),\r
  83   fTotalNumberOfAddedNtupleEntries(0)\r
  84 {\r
  85   // Default constructor\r
  86   // Common I/O in slot 0\r
  87   DefineInput (0, TChain::Class());\r
  88   DefineOutput(0, TTree::Class());\r
  89         \r
  90   // Your private output\r
  91   DefineOutput(1, TList::Class());\r
  92 }\r
  93 \r
  94 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
  95   AliAnalysisTaskSE(name),\r
  96   fV0Reader(NULL),\r
  97   fStack(NULL),\r
  98   fESDEvent(NULL),      \r
  99   fOutputContainer(0x0),\r
 100   fHistograms(NULL),\r
 101   fDoMCTruth(kFALSE),\r
 102   fMCAllGammas(),\r
 103   fMCPi0s(),\r
 104   fMCEtas(),\r
 105   fMCGammaChic(),\r
 106   fKFReconstructedGammas(),\r
 107   fIsTrueReconstructedGammas(),\r
 108   fElectronv1(),\r
 109   fElectronv2(),\r
 110   fCurrentEventPosElectron(),\r
 111   fPreviousEventPosElectron(),\r
 112   fCurrentEventNegElectron(),\r
 113   fPreviousEventNegElectron(),\r
 114   fKFReconstructedGammasCut(),  \r
 115   fPreviousEventTLVNegElectron(),\r
 116   fPreviousEventTLVPosElectron(),\r
 117   fElectronMass(-1),\r
 118   fGammaMass(-1),\r
 119   fPi0Mass(-1),\r
 120   fEtaMass(-1),\r
 121   fGammaWidth(-1),\r
 122   fPi0Width(-1),\r
 123   fEtaWidth(-1),\r
 124   fMinOpeningAngleGhostCut(0.),\r
 125   fCalculateBackground(kFALSE),\r
 126   fWriteNtuple(kFALSE),\r
 127   fGammaNtuple(NULL),\r
 128   fNeutralMesonNtuple(NULL),\r
 129   fTotalNumberOfAddedNtupleEntries(0)\r
 130 {\r
 131   // Common I/O in slot 0\r
 132   DefineInput (0, TChain::Class());\r
 133   DefineOutput(0, TTree::Class());\r
 134         \r
 135   // Your private output\r
 136   DefineOutput(1, TList::Class());\r
 137 }\r
 138 \r
 139 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
 140 {\r
 141   // Remove all pointers\r
 142         \r
 143   if(fOutputContainer){\r
 144     fOutputContainer->Clear() ; \r
 145     delete fOutputContainer ;\r
 146   }\r
 147   if(fHistograms){\r
 148     delete fHistograms;\r
 149   }\r
 150   if(fV0Reader){\r
 151     delete fV0Reader;\r
 152   }\r
 153 }\r
 154 \r
 155 \r
 156 void AliAnalysisTaskGammaConversion::Init()\r
 157 {\r
 158   // Initialization\r
 159   AliLog::SetGlobalLogLevel(AliLog::kError);\r
 160 }\r
 161 \r
 162 \r
 163 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
 164 {\r
 165   // Execute analysis for current event\r
 166         \r
 167   ConnectInputData("");\r
 168         \r
 169   //clear vectors\r
 170   fMCAllGammas.clear();\r
 171   fMCPi0s.clear();\r
 172   fMCEtas.clear();\r
 173   fMCGammaChic.clear();\r
 174         \r
 175   fKFReconstructedGammas.clear();\r
 176   fIsTrueReconstructedGammas.clear();\r
 177   fElectronv1.clear();\r
 178   fElectronv2.clear();\r
 179   fCurrentEventPosElectron.clear();\r
 180   fCurrentEventNegElectron.clear();     \r
 181   fKFReconstructedGammasCut.clear(); \r
 182         \r
 183   //Clear the data in the v0Reader\r
 184   fV0Reader->UpdateEventByEventData();\r
 185 \r
 186   \r
 187   // Process the MC information\r
 188   if(fDoMCTruth){\r
 189     ProcessMCData();\r
 190   }\r
 191   \r
 192   //Process the v0 information with no cuts\r
 193   ProcessV0sNoCut();\r
 194   \r
 195   // Process the v0 information\r
 196   ProcessV0s();\r
 197   \r
 198   //calculate background if flag is set\r
 199   if(fCalculateBackground){\r
 200     CalculateBackground();\r
 201   }\r
 202   \r
 203   // Process reconstructed gammas\r
 204   ProcessGammasForNeutralMesonAnalysis();\r
 205 \r
 206   CheckV0Efficiency();\r
 207   \r
 208   //Process reconstructed gammas electrons for Chi_c Analysis\r
 209   ProcessGammaElectronsForChicAnalysis();\r
 210 \r
 211   PostData(1, fOutputContainer);\r
 212         \r
 213 }\r
 214 \r
 215 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
 216   // see header file for documentation\r
 217         \r
 218   if(fV0Reader == NULL){\r
 219     // Write warning here cuts and so on are default if this ever happens\r
 220   }\r
 221   fV0Reader->Initialize();\r
 222 }\r
 223 \r
 224 \r
 225 \r
 226 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
 227   // see header file for documentation\r
 228         \r
 229   fStack = fV0Reader->GetMCStack();\r
 230 \r
 231   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
 232     return; // aborts if the primary vertex does not have contributors.\r
 233   }\r
 234 \r
 235   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
 236     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
 237                 \r
 238     if (!particle) {\r
 239       //print warning here\r
 240       continue;\r
 241     }\r
 242 \r
 243     ///////////////////////Begin Chic Analysis/////////////////////////////\r
 244 \r
 245 \r
 246     if(particle->GetPdgCode() == 443){//Is JPsi\r
 247 \r
 248       if(particle->GetNDaughters()==2){\r
 249         if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&\r
 250            TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){\r
 251           TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
 252           TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
 253           if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)\r
 254             fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance\r
 255 \r
 256           if( TMath::Abs(daug0->Eta()) < 0.9){\r
 257             if(daug0->GetPdgCode() == -11)\r
 258               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance\r
 259             else\r
 260               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance\r
 261 \r
 262           }\r
 263           if(TMath::Abs(daug1->Eta()) < 0.9){\r
 264             if(daug1->GetPdgCode() == -11)\r
 265               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance\r
 266             else\r
 267               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance\r
 268           }\r
 269         }\r
 270       }\r
 271     }\r
 272     //              const int CHI_C0   = 10441;\r
 273     //              const int CHI_C1   = 20443;\r
 274     //              const int CHI_C2   = 445\r
 275     if(particle->GetPdgCode() == 22){//gamma from JPsi\r
 276       if(particle->GetMother(0) > -1){\r
 277         if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||\r
 278            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||\r
 279            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){\r
 280           if(TMath::Abs(particle->Eta()) < 1.2)\r
 281             fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance\r
 282         }\r
 283       }\r
 284     }\r
 285     if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){\r
 286       if( particle->GetNDaughters() == 2){\r
 287         TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
 288         TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
 289 \r
 290         if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){\r
 291           if( daug0->GetPdgCode() == 443){\r
 292             TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());\r
 293             TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());\r
 294             if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
 295               fHistograms->FillTable("Table_Electrons",18);\r
 296 \r
 297           }//if\r
 298           else if (daug1->GetPdgCode() == 443){\r
 299             TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());\r
 300             TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());\r
 301             if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
 302               fHistograms->FillTable("Table_Electrons",18);\r
 303           }//else if\r
 304         }//gamma o Jpsi\r
 305       }//GetNDaughters\r
 306     }\r
 307 \r
 308 \r
 309     /////////////////////End Chic Analysis////////////////////////////\r
 310 \r
 311 \r
 312     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;\r
 313                                         \r
 314     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point\r
 315                 \r
 316     Double_t tmpPhi=particle->Phi();\r
 317                 \r
 318     if(particle->Phi()> TMath::Pi()){\r
 319       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
 320     }\r
 321                 \r
 322     Double_t rapidity;\r
 323     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){\r
 324       rapidity=0;\r
 325     }\r
 326     else{\r
 327       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));\r
 328     }   \r
 329                 \r
 330     //process the gammas\r
 331     if (particle->GetPdgCode() == 22){\r
 332                         \r
 333       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){\r
 334         continue; // no photon as mothers!\r
 335       }\r
 336 \r
 337       if(particle->GetMother(0) >= fStack->GetNprimary()){\r
 338         continue; // the gamma has a mother, and it is not a primary particle\r
 339       }\r
 340 \r
 341       fMCAllGammas.push_back(particle);\r
 342                         \r
 343       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());\r
 344       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());\r
 345       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());\r
 346       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);\r
 347       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);\r
 348                         \r
 349                         \r
 350       if(particle->GetMother(0) < 0){   // direct gamma\r
 351         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());\r
 352         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());\r
 353         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());\r
 354         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);\r
 355         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);                                \r
 356       }\r
 357                         \r
 358                         \r
 359       // looking for conversion (electron + positron from pairbuilding (= 5) )\r
 360       TParticle* ePos = NULL;\r
 361       TParticle* eNeg = NULL;\r
 362                         \r
 363       if(particle->GetNDaughters() >= 2){\r
 364         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
 365           TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
 366           if(tmpDaughter->GetUniqueID() == 5){\r
 367             if(tmpDaughter->GetPdgCode() == 11){\r
 368               eNeg = tmpDaughter;\r
 369             }\r
 370             else if(tmpDaughter->GetPdgCode() == -11){\r
 371               ePos = tmpDaughter;\r
 372             }\r
 373           }\r
 374         }\r
 375       }\r
 376                         \r
 377                         \r
 378       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production\r
 379         continue;\r
 380       }\r
 381                         \r
 382                         \r
 383       Double_t ePosPhi = ePos->Phi();\r
 384       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());\r
 385                         \r
 386       Double_t eNegPhi = eNeg->Phi();\r
 387       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());\r
 388                         \r
 389                         \r
 390       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){\r
 391         continue; // no reconstruction below the Pt cut\r
 392       }\r
 393                                         \r
 394       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){\r
 395         continue;\r
 396       } \r
 397                                 \r
 398       if(ePos->R()>fV0Reader->GetMaxRCut()){\r
 399         continue; // cuts on distance from collision point\r
 400       }\r
 401       \r
 402       \r
 403       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){\r
 404         continue;               // line cut to exclude regions where we do not reconstruct\r
 405       }         \r
 406                 \r
 407       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());\r
 408       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());\r
 409       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());\r
 410       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);\r
 411       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);\r
 412       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());\r
 413                         \r
 414       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());\r
 415       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());\r
 416       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());\r
 417       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);\r
 418                         \r
 419       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());\r
 420       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());\r
 421       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());\r
 422       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);\r
 423                         \r
 424                         \r
 425                         \r
 426       //cout << "filled histos for converted gamma, ePos, eNeg" << endl;\r
 427                         \r
 428       // begin Mapping \r
 429       Int_t rBin    = fHistograms->GetRBin(ePos->R());\r
 430       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());\r
 431                         \r
 432       TString nameMCMappingPhiR="";\r
 433       nameMCMappingPhiR.Form("MC_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 434       fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());\r
 435                         \r
 436       TString nameMCMappingPhi="";\r
 437       nameMCMappingPhi.Form("MC_Conversion_Mapping-Phi%02d",phiBin);\r
 438       fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
 439                         \r
 440       TString nameMCMappingR="";\r
 441       nameMCMappingR.Form("MC_Conversion_Mapping-R%02d",rBin);\r
 442       fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
 443                         \r
 444       TString nameMCMappingPhiInR="";\r
 445       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_R-%02d",rBin);\r
 446       fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
 447       //end mapping\r
 448                         \r
 449       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());\r
 450       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());\r
 451       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());\r
 452       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));\r
 453                         \r
 454       //cout << "mapping is done" << endl;\r
 455                         \r
 456                         \r
 457       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted\r
 458         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());\r
 459         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());\r
 460         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());\r
 461         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);\r
 462         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);\r
 463                                 \r
 464       } // end direct gamma\r
 465       else{   // mother exits \r
 466         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
 467             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
 468             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
 469             ){ \r
 470           fMCGammaChic.push_back(particle);\r
 471         }\r
 472       }  // end if mother exits\r
 473     } // end if particle is a photon\r
 474                 \r
 475     if(particle->GetNDaughters() == 2){\r
 476                         \r
 477       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
 478       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
 479                         \r
 480       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters\r
 481                         \r
 482                         \r
 483                         \r
 484       // check for conversions now -> have to pass eta and line cut!\r
 485       Bool_t daughter0Electron = kFALSE;\r
 486       Bool_t daughter0Positron = kFALSE;\r
 487       Bool_t daughter1Electron = kFALSE;\r
 488       Bool_t daughter1Positron = kFALSE;\r
 489                         \r
 490                         \r
 491                         \r
 492       if(daughter0->GetNDaughters() >= 2){\r
 493         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){\r
 494           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
 495           if(tmpDaughter->GetUniqueID() == 5){\r
 496             if(tmpDaughter->GetPdgCode() == 11){\r
 497               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 498                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 499                   daughter0Electron = kTRUE;\r
 500                 }\r
 501                                                                 \r
 502               }\r
 503             }\r
 504             else if(tmpDaughter->GetPdgCode() == -11){\r
 505               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 506                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 507                   daughter0Positron = kTRUE;\r
 508                 }\r
 509                                                                 \r
 510               }\r
 511                                                         \r
 512             }\r
 513           }\r
 514         }\r
 515       }\r
 516                         \r
 517                         \r
 518         \r
 519       if(daughter1->GetNDaughters() >= 2){\r
 520         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){\r
 521           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
 522           if(tmpDaughter->GetUniqueID() == 5){\r
 523             if(tmpDaughter->GetPdgCode() == 11){\r
 524               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 525                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 526                   daughter1Electron = kTRUE;\r
 527                 }\r
 528                                                                 \r
 529               }\r
 530             }\r
 531             else if(tmpDaughter->GetPdgCode() == -11){\r
 532               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 533                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 534                   daughter1Positron = kTRUE;\r
 535                 }\r
 536                                                                 \r
 537               }\r
 538                                                         \r
 539             }\r
 540           }\r
 541         }\r
 542       }\r
 543                         \r
 544                                                                                                 \r
 545                         \r
 546                         \r
 547       if(particle->GetPdgCode()==111){     //Pi0\r
 548         if( iTracks >= fStack->GetNprimary()){\r
 549           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
 550           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);\r
 551           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
 552           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
 553           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
 554           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
 555           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
 556           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 557           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
 558                                         \r
 559           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 560             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 561             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 562             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 563               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 564               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 565             }\r
 566           }\r
 567         }\r
 568         else{\r
 569           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());    \r
 570           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);\r
 571           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
 572           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
 573           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
 574           fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
 575           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
 576           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 577           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 578                                         \r
 579           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 580             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 581             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 582             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 583               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 584               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 585             }\r
 586           }\r
 587         }\r
 588       }\r
 589                         \r
 590       if(particle->GetPdgCode()==221){   //Eta\r
 591         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
 592         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);\r
 593         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
 594         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
 595         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
 596         fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
 597         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
 598         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 599         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 600                                 \r
 601         if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 602           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 603           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 604           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 605             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 606             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 607           }\r
 608                                         \r
 609         }\r
 610                                 \r
 611       }\r
 612                         \r
 613       // all motherparticles with 2 gammas as daughters\r
 614       fHistograms->FillHistogram("MC_Mother_R", particle->R());\r
 615       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());\r
 616       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());\r
 617       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());\r
 618       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 619       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());\r
 620       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());\r
 621       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());\r
 622       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);\r
 623       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);\r
 624       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());                 \r
 625       if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 626         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 627         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 628         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());                      \r
 629         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 630           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 631           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 632           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());                  \r
 633 \r
 634         }\r
 635                                 \r
 636                                 \r
 637       }\r
 638                         \r
 639       //cout << "mother histos are filled" << endl;\r
 640                         \r
 641     } // end if(particle->GetNDaughters() == 2)\r
 642                 \r
 643   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
 644         \r
 645   //cout << "right before the end of processMCdata" << endl;\r
 646         \r
 647 } // end ProcessMCData\r
 648 \r
 649 \r
 650 \r
 651 void AliAnalysisTaskGammaConversion::FillNtuple(){\r
 652   //Fills the ntuple with the different values\r
 653 \r
 654   if(fGammaNtuple == NULL){\r
 655     return;\r
 656   }\r
 657   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
 658   for(Int_t i=0;i<numberOfV0s;i++){\r
 659     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
 660     AliESDv0 * cV0 = fV0Reader->GetV0(i);\r
 661     Double_t negPID=0;\r
 662     Double_t posPID=0;\r
 663     fV0Reader->GetPIDProbability(negPID,posPID);\r
 664     values[0]=cV0->GetOnFlyStatus();\r
 665     values[1]=fV0Reader->CheckForPrimaryVertex();\r
 666     values[2]=negPID;\r
 667     values[3]=posPID;\r
 668     values[4]=fV0Reader->GetX();\r
 669     values[5]=fV0Reader->GetY();\r
 670     values[6]=fV0Reader->GetZ();\r
 671     values[7]=fV0Reader->GetXYRadius();\r
 672     values[8]=fV0Reader->GetMotherCandidateNDF();\r
 673     values[9]=fV0Reader->GetMotherCandidateChi2();\r
 674     values[10]=fV0Reader->GetMotherCandidateEnergy();\r
 675     values[11]=fV0Reader->GetMotherCandidateEta();\r
 676     values[12]=fV0Reader->GetMotherCandidatePt();\r
 677     values[13]=fV0Reader->GetMotherCandidateMass();\r
 678     values[14]=fV0Reader->GetMotherCandidateWidth();\r
 679     //    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
 680     values[16]=fV0Reader->GetOpeningAngle();\r
 681     values[17]=fV0Reader->GetNegativeTrackEnergy();\r
 682     values[18]=fV0Reader->GetNegativeTrackPt();\r
 683     values[19]=fV0Reader->GetNegativeTrackEta();\r
 684     values[20]=fV0Reader->GetNegativeTrackPhi();\r
 685     values[21]=fV0Reader->GetPositiveTrackEnergy();\r
 686     values[22]=fV0Reader->GetPositiveTrackPt();\r
 687     values[23]=fV0Reader->GetPositiveTrackEta();\r
 688     values[24]=fV0Reader->GetPositiveTrackPhi();\r
 689     values[25]=fV0Reader->HasSameMCMother();\r
 690     if(values[25] != 0){\r
 691       values[26]=fV0Reader->GetMotherMCParticlePDGCode();\r
 692       values[15]=fV0Reader->GetMotherMCParticle()->Pt();\r
 693     }\r
 694     fTotalNumberOfAddedNtupleEntries++;\r
 695     fGammaNtuple->Fill(values);\r
 696   }\r
 697   fV0Reader->ResetV0IndexNumber();\r
 698         \r
 699 }\r
 700 \r
 701 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){\r
 702   // Process all the V0's without applying any cuts to it\r
 703 \r
 704   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
 705   for(Int_t i=0;i<numberOfV0s;i++){\r
 706     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);\r
 707 \r
 708     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
 709       return;\r
 710     }\r
 711     \r
 712     if(fDoMCTruth){\r
 713       \r
 714       if(fV0Reader->HasSameMCMother() == kFALSE){\r
 715         continue;\r
 716       }\r
 717                 \r
 718       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
 719       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
 720 \r
 721       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
 722         continue;\r
 723       }\r
 724       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
 725         continue;\r
 726       }\r
 727         \r
 728       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
 729       \r
 730         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 731         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 732         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                               \r
 733         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 734         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 735         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 736         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 737         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 738         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 739         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 740         \r
 741         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 742         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());\r
 743         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 744         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 745         \r
 746         /*\r
 747           ESD_NoCutConvGamma_Pt\r
 748           ESD_NoCutConvGamma_Energy\r
 749           ESD_NoCutConvGamma_Eta\r
 750           ESD_NoCutConvGamma_Phi\r
 751           ESD_NoCutConvGamma_Mass\r
 752           ESD_NoCutConvGamma_Width\r
 753           ESD_NoCutConvGamma_Chi2\r
 754           ESD_NoCutConvGamma_NDF\r
 755           ESD_NoCutConvGamma_PtvsEta\r
 756           ESD_NoCutConversion_XY\r
 757           ESD_NoCutConversion_R\r
 758           ESD_NoCutConversion_ZR\r
 759           ESD_NoCutConversion_OpeningAngle\r
 760         */\r
 761       }\r
 762     }\r
 763   }\r
 764   fV0Reader->ResetV0IndexNumber();\r
 765 }\r
 766 \r
 767 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
 768   // see header file for documentation\r
 769         \r
 770   if(fWriteNtuple == kTRUE){\r
 771     FillNtuple();\r
 772   }\r
 773         \r
 774   Int_t nSurvivingV0s=0;\r
 775   while(fV0Reader->NextV0()){\r
 776     nSurvivingV0s++;\r
 777                 \r
 778                 \r
 779     //-------------------------- filling v0 information -------------------------------------\r
 780     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());\r
 781     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 782     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 783     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
 784                 \r
 785     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
 786     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
 787     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
 788     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
 789                 \r
 790     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
 791     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
 792     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
 793     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
 794                 \r
 795     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 796     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 797     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());\r
 798     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 799     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 800     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 801     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 802     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 803     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 804     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 805                 \r
 806                 \r
 807     // begin mapping\r
 808     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
 809     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
 810     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
 811                 \r
 812     TString nameESDMappingPhiR="";\r
 813     nameESDMappingPhiR.Form("ESD_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 814     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
 815                 \r
 816     TString nameESDMappingPhi="";\r
 817     nameESDMappingPhi.Form("ESD_Conversion_Mapping-Phi%02d",phiBin);\r
 818     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
 819                 \r
 820     TString nameESDMappingR="";\r
 821     nameESDMappingR.Form("ESD_Conversion_Mapping-R%02d",rBin);\r
 822     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
 823                 \r
 824     TString nameESDMappingPhiInR="";\r
 825     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_R-%02d",rBin);\r
 826     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
 827     // end mapping\r
 828                 \r
 829     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
 830     fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());\r
 831     fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());\r
 832 \r
 833                 \r
 834     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
 835     if(fDoMCTruth){\r
 836                         \r
 837       if(fV0Reader->HasSameMCMother() == kFALSE){\r
 838         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 839         continue;\r
 840       }\r
 841       \r
 842                 \r
 843       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
 844       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
 845 \r
 846       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
 847         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 848         continue;\r
 849       }\r
 850       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
 851         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 852         continue;\r
 853       }\r
 854         \r
 855       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
 856         fIsTrueReconstructedGammas.push_back(kTRUE);\r
 857                                 \r
 858         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 859         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 860         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                                \r
 861         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 862         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 863         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 864         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 865         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 866         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 867         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 868         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());\r
 869         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());\r
 870         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
 871         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
 872         \r
 873         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 874         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
 875         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 876         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 877 \r
 878         \r
 879         /*\r
 880           fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 881           fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
 882           fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 883           fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 884         */\r
 885 \r
 886 \r
 887 \r
 888         //resolution\r
 889         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
 890         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
 891         Double_t resdPt = 0;\r
 892         if(mcpt != 0){\r
 893           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
 894         }\r
 895                                 \r
 896         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
 897         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
 898         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
 899                                 \r
 900         Double_t resdZ = 0;\r
 901         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
 902           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
 903         }\r
 904                                 \r
 905         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
 906         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
 907         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
 908                                 \r
 909         Double_t resdR = 0;\r
 910         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
 911           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
 912         }\r
 913 \r
 914         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
 915         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
 916         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
 917         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
 918       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)\r
 919       else{\r
 920         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 921       }\r
 922     }//if(fDoMCTruth)\r
 923   }//while(fV0Reader->NextV0)\r
 924   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);\r
 925   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
 926         \r
 927   //cout << "nearly at the end of doMCTruth" << endl;\r
 928         \r
 929 }\r
 930 \r
 931 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
 932   // see header file for documentation\r
 933         \r
 934   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
 935     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
 936                         \r
 937       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
 938       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
 939       \r
 940       if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
 941         continue;\r
 942       }\r
 943       if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
 944         continue;\r
 945       }\r
 946 \r
 947       /*\r
 948         if(fIsTrueReconstructedGammas[firstGammaIndex] == kFALSE || fIsTrueReconstructedGammas[secondGammaIndex] == kFALSE){\r
 949         continue;\r
 950         }\r
 951       */\r
 952 \r
 953       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
 954                         \r
 955       Double_t massTwoGammaCandidate = 0.;\r
 956       Double_t widthTwoGammaCandidate = 0.;\r
 957       Double_t chi2TwoGammaCandidate =10000.;   \r
 958       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
 959       if(twoGammaCandidate->GetNDF()>0){\r
 960         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
 961                                 \r
 962         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){                                      \r
 963                                         \r
 964           TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());\r
 965           TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());\r
 966                                         \r
 967           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);                                 \r
 968           Double_t rapidity;\r
 969           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){\r
 970             rapidity=0;\r
 971           }\r
 972           else{\r
 973             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));\r
 974           }\r
 975                                         \r
 976           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut) continue;   // minimum opening angle to avoid using ghosttracks\r
 977 \r
 978           fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
 979           fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());\r
 980           fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());\r
 981           fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());\r
 982           fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);                                     \r
 983           fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());\r
 984           fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);\r
 985           fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!\r
 986           fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());\r
 987           fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());\r
 988           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());\r
 989           fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);\r
 990         }\r
 991       }\r
 992       delete twoGammaCandidate;\r
 993                         \r
 994       //cout << "nearly at the end of processgamma for neutral meson ..." << endl;\r
 995                         \r
 996                         \r
 997     }\r
 998   }\r
 999 }\r
1000 \r
1001 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
1002   // see header file for documentation\r
1003         \r
1004   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
1005   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
1006   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
1007     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
1008     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
1009       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
1010 \r
1011       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
1012                         \r
1013       Double_t massBG =0.;\r
1014       Double_t widthBG = 0.;\r
1015       Double_t chi2BG =10000.;  \r
1016       backgroundCandidate->GetMass(massBG,widthBG);\r
1017       if(backgroundCandidate->GetNDF()>0){\r
1018         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
1019         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
1020                                         \r
1021           TVector3 MomentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());\r
1022           TVector3 SpaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());\r
1023                                         \r
1024           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
1025 \r
1026           Double_t rapidity;\r
1027           if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;\r
1028           else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));\r
1029 \r
1030                                         \r
1031                                         \r
1032                                         \r
1033           if(openingAngleBG < fMinOpeningAngleGhostCut ) continue;   // minimum opening angle to avoid using ghosttracks\r
1034                         \r
1035                                         \r
1036           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
1037           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
1038           fHistograms->FillHistogram("ESD_Background_Pt",  MomentumVectorbackgroundCandidate.Pt());\r
1039           fHistograms->FillHistogram("ESD_Background_Eta", MomentumVectorbackgroundCandidate.Eta());\r
1040           fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);\r
1041           fHistograms->FillHistogram("ESD_Background_Phi", SpaceVectorbackgroundCandidate.Phi());\r
1042           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
1043           fHistograms->FillHistogram("ESD_Background_R", SpaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!\r
1044           fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());\r
1045           fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());\r
1046           fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());\r
1047           fHistograms->FillHistogram("ESD_Background_InvMass",massBG);\r
1048         }\r
1049       }\r
1050       delete backgroundCandidate;   \r
1051       //cout << "nearly at the end of background" << endl;\r
1052                         \r
1053     }\r
1054   }\r
1055 }\r
1056 \r
1057 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
1058 {\r
1059   // Terminate analysis\r
1060   //\r
1061   AliDebug(1,"Do nothing in Terminate");\r
1062 }\r
1063 \r
1064 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
1065 {\r
1066   // Create the output container\r
1067   if(fOutputContainer != NULL){\r
1068     delete fOutputContainer;\r
1069     fOutputContainer = NULL;\r
1070   }\r
1071   if(fOutputContainer == NULL){\r
1072     fOutputContainer = new TList();\r
1073   }\r
1074         \r
1075   //Adding the histograms to the output container\r
1076   fHistograms->GetOutputContainer(fOutputContainer);\r
1077         \r
1078         \r
1079   if(fWriteNtuple){\r
1080     if(fGammaNtuple == NULL){\r
1081       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
1082     }\r
1083     if(fNeutralMesonNtuple == NULL){\r
1084       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");\r
1085     }\r
1086     TList * ntupleTList = new TList();\r
1087     ntupleTList->SetName("Ntuple");\r
1088     ntupleTList->Add((TNtuple*)fGammaNtuple);\r
1089     fOutputContainer->Add(ntupleTList);\r
1090   }\r
1091         \r
1092   fOutputContainer->SetName(GetName());\r
1093 }\r
1094 \r
1095 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{\r
1096   //helper function\r
1097   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
1098   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
1099   return v3D0.Angle(v3D1);\r
1100 }\r
1101 \r
1102 void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){\r
1103   // see header file for documentation\r
1104 \r
1105   vector<Int_t> indexOfGammaParticle;\r
1106 \r
1107   fStack = fV0Reader->GetMCStack();\r
1108 \r
1109   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
1110     return; // aborts if the primary vertex does not have contributors.\r
1111   }\r
1112 \r
1113   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {\r
1114     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
1115     if(particle->GetPdgCode()==22){     //Gamma\r
1116       if(particle->GetNDaughters() >= 2){\r
1117         TParticle* electron=NULL;\r
1118         TParticle* positron=NULL; \r
1119         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
1120           TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
1121           if(tmpDaughter->GetUniqueID() == 5){\r
1122             if(tmpDaughter->GetPdgCode() == 11){\r
1123               electron = tmpDaughter;\r
1124             }\r
1125             else if(tmpDaughter->GetPdgCode() == -11){\r
1126               positron = tmpDaughter;\r
1127             }\r
1128           }\r
1129         }\r
1130         if(electron!=NULL && positron!=0){\r
1131           if(electron->R()<160){\r
1132             indexOfGammaParticle.push_back(iTracks);\r
1133           }\r
1134         }\r
1135       }\r
1136     }\r
1137   }\r
1138 \r
1139   Int_t nFoundGammas=0;\r
1140   Int_t nNotFoundGammas=0;\r
1141 \r
1142   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
1143   for(Int_t i=0;i<numberOfV0s;i++){\r
1144     fV0Reader->GetV0(i);\r
1145     \r
1146     if(fV0Reader->HasSameMCMother() == kFALSE){\r
1147       continue;\r
1148     }\r
1149     \r
1150     TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
1151     TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
1152     \r
1153     if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
1154       continue;\r
1155     }\r
1156     if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
1157       continue;\r
1158     }\r
1159     \r
1160     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
1161       //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();\r
1162       for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){\r
1163         if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){\r
1164           nFoundGammas++;\r
1165         }\r
1166         else{\r
1167           nNotFoundGammas++;\r
1168         }\r
1169       }\r
1170     }\r
1171   }\r
1172   //  cout<<"Found: "<<nFoundGammas<<"  of: "<<indexOfGammaParticle.size()<<endl;\r
1173 }\r
1174 \r
1175 \r
1176 void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){\r
1177   // see header file for documantation\r
1178 \r
1179   fESDEvent = fV0Reader->GetESDEvent();\r
1180 \r
1181 \r
1182   vector <AliESDtrack*> vESDeNegTemp(0);\r
1183   vector <AliESDtrack*> vESDePosTemp(0);\r
1184   vector <AliESDtrack*> vESDxNegTemp(0);\r
1185   vector <AliESDtrack*> vESDxPosTemp(0);\r
1186   vector <AliESDtrack*> vESDeNegNoJPsi(0);\r
1187   vector <AliESDtrack*> vESDePosNoJPsi(0); \r
1188 \r
1189 \r
1190 \r
1191   fHistograms->FillTable("Table_Electrons",0);//Count number of Events\r
1192 \r
1193   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){\r
1194     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);\r
1195 \r
1196     if(!curTrack){\r
1197       //print warning here\r
1198       continue;\r
1199     }\r
1200 \r
1201     double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;\r
1202     double r[3];curTrack->GetConstrainedXYZ(r);\r
1203 \r
1204     TVector3 rXYZ(r);\r
1205 \r
1206     fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks\r
1207 \r
1208     Bool_t flagKink       =  kTRUE;\r
1209     Bool_t flagTPCrefit   =  kTRUE;\r
1210     Bool_t flagTRDrefit   =  kTRUE;\r
1211     Bool_t flagITSrefit   =  kTRUE;\r
1212     Bool_t flagTRDout     =  kTRUE;\r
1213     Bool_t flagVertex     =  kTRUE;\r
1214 \r
1215 \r
1216     //Cuts ---------------------------------------------------------------\r
1217 \r
1218     if(curTrack->GetKinkIndex(0) > 0){\r
1219       fHistograms->FillHistogram("Table_Electrons",5);//Count kink\r
1220       flagKink = kFALSE;\r
1221     }\r
1222 \r
1223     ULong_t trkStatus = curTrack->GetStatus();\r
1224 \r
1225     ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);\r
1226 \r
1227     if(!tpcRefit){\r
1228       fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit\r
1229       flagTPCrefit = kFALSE;\r
1230     }\r
1231 \r
1232     ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);\r
1233     if(!itsRefit){\r
1234       fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit\r
1235       flagITSrefit = kFALSE;\r
1236     }\r
1237 \r
1238     ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);\r
1239 \r
1240     if(!trdRefit){\r
1241       fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit\r
1242       flagTRDrefit = kFALSE;\r
1243     }\r
1244 \r
1245     ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);\r
1246 \r
1247     if(!trdOut) {\r
1248       fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout\r
1249       flagTRDout = kFALSE;\r
1250     }\r
1251 \r
1252     double nsigmaToVxt = GetSigmaToVertex(curTrack);\r
1253 \r
1254     if(nsigmaToVxt > 3){\r
1255       fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3\r
1256       flagVertex = kFALSE;\r
1257     }\r
1258 \r
1259     if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;\r
1260     fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts\r
1261 \r
1262 \r
1263     Stat_t pid, weight;\r
1264     GetPID(curTrack, pid, weight);\r
1265 \r
1266     if(pid!=0){\r
1267       fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0\r
1268     }\r
1269 \r
1270     if(pid == 0){\r
1271       fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0\r
1272     }\r
1273 \r
1274 \r
1275 \r
1276 \r
1277     Int_t labelMC = TMath::Abs(curTrack->GetLabel());\r
1278     TParticle* curParticle = fStack->Particle(labelMC);\r
1279 \r
1280 \r
1281 \r
1282 \r
1283     TLorentzVector curElec;\r
1284     curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);\r
1285 \r
1286 \r
1287 \r
1288 \r
1289     if(curTrack->GetSign() > 0){\r
1290 \r
1291       vESDxPosTemp.push_back(curTrack);\r
1292 \r
1293       if( pid == 0){\r
1294 \r
1295         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
1296         fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());\r
1297         fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
1298         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
1299         fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
1300         vESDePosTemp.push_back(curTrack);\r
1301 \r
1302 \r
1303 \r
1304       }\r
1305 \r
1306     }\r
1307     else {\r
1308       vESDxNegTemp.push_back(curTrack);\r
1309 \r
1310       if( pid == 0){\r
1311 \r
1312         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
1313         fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());\r
1314         fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
1315         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
1316         fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
1317         vESDeNegTemp.push_back(curTrack);\r
1318 \r
1319 \r
1320 \r
1321 \r
1322       }\r
1323     }\r
1324 \r
1325   }\r
1326 \r
1327 \r
1328   Bool_t ePosJPsi = kFALSE;\r
1329   Bool_t eNegJPsi = kFALSE;             \r
1330   Bool_t ePosPi0  = kFALSE;\r
1331   Bool_t eNegPi0  = kFALSE;\r
1332         \r
1333   UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;\r
1334  \r
1335   for(UInt_t iNeg=0; iNeg < vESDeNegTemp.size(); iNeg++){\r
1336     if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetPdgCode() == 11)\r
1337       if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0) > -1){\r
1338         Int_t labelMother = fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0);\r
1339         TParticle* partMother = fStack ->Particle(labelMother);\r
1340         if (partMother->GetPdgCode() == 111){\r
1341           ieNegPi0 = iNeg;\r
1342           eNegPi0 = kTRUE;\r
1343         }\r
1344         if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
1345           fHistograms->FillTable("Table_Electrons",14);\r
1346           ieNegJPsi = iNeg;\r
1347           eNegJPsi = kTRUE;\r
1348         }\r
1349         else{   \r
1350           vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);\r
1351           //            cout<<"ESD No Positivo JPsi "<<endl;\r
1352         }\r
1353 \r
1354       }\r
1355   }     \r
1356 \r
1357   for(UInt_t iPos=0; iPos < vESDePosTemp.size(); iPos++){\r
1358     if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetPdgCode() == -11)\r
1359       if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0) > -1){\r
1360         Int_t labelMother = fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0);\r
1361         TParticle* partMother = fStack ->Particle(labelMother);\r
1362         if (partMother->GetPdgCode() == 111){\r
1363           iePosPi0 = iPos;\r
1364           ePosPi0 = kTRUE;\r
1365         }\r
1366         if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
1367           fHistograms->FillTable("Table_Electrons",15);\r
1368           iePosJPsi = iPos;\r
1369           ePosJPsi = kTRUE;\r
1370         }\r
1371         else{\r
1372           vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);\r
1373           //            cout<<"ESD No Negativo JPsi "<<endl;\r
1374         }\r
1375 \r
1376       }\r
1377   }\r
1378         \r
1379   if( eNegJPsi && ePosJPsi ){\r
1380     TVector3 tempeNegV,tempePosV;\r
1381     tempeNegV.SetXYZ(vESDeNegTemp[ieNegJPsi]->Px(),vESDeNegTemp[ieNegJPsi]->Py(),vESDeNegTemp[ieNegJPsi]->Pz());                        \r
1382     tempePosV.SetXYZ(vESDePosTemp[iePosJPsi]->Px(),vESDePosTemp[iePosJPsi]->Py(),vESDePosTemp[iePosJPsi]->Pz());\r
1383     fHistograms->FillTable("Table_Electrons",16);\r
1384     fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));       \r
1385     fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegJPsi]->GetLabel())),\r
1386                                                                               fStack->Particle(TMath::Abs(vESDePosTemp[iePosJPsi]->GetLabel()))));      \r
1387   }\r
1388         \r
1389   if( eNegPi0 && ePosPi0 ){\r
1390     TVector3 tempeNegV,tempePosV;\r
1391     tempeNegV.SetXYZ(vESDeNegTemp[ieNegPi0]->Px(),vESDeNegTemp[ieNegPi0]->Py(),vESDeNegTemp[ieNegPi0]->Pz());\r
1392     tempePosV.SetXYZ(vESDePosTemp[iePosPi0]->Px(),vESDePosTemp[iePosPi0]->Py(),vESDePosTemp[iePosPi0]->Pz());\r
1393     fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));\r
1394     fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegPi0]->GetLabel())),\r
1395                                                                              fStack->Particle(TMath::Abs(vESDePosTemp[iePosPi0]->GetLabel()))));   \r
1396   }\r
1397          \r
1398 \r
1399   FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));\r
1400 \r
1401   CleanWithAngleCuts(vESDeNegTemp,vESDePosTemp,fKFReconstructedGammas);\r
1402         \r
1403   vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);\r
1404   vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);\r
1405 \r
1406 \r
1407   FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);\r
1408 \r
1409  \r
1410 \r
1411 \r
1412   //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);\r
1413 \r
1414 \r
1415   FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);\r
1416   FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));\r
1417 \r
1418        \r
1419 \r
1420   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",\r
1421                            fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);\r
1422 \r
1423   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",\r
1424                            fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);\r
1425 \r
1426   //BackGround\r
1427 \r
1428   //Like Sign e+e-\r
1429   ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);\r
1430   ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);\r
1431   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);\r
1432   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);\r
1433 \r
1434   //        Like Sign e+e- no JPsi\r
1435   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));\r
1436   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));\r
1437 \r
1438   //Mixed Event\r
1439 \r
1440   if( fCurrentEventPosElectron.size() > 0 && fCurrentEventNegElectron.size() > 0 && fKFReconstructedGammasCut.size() > 0 ){\r
1441     FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",\r
1442                              fKFReconstructedGammasCut,fPreviousEventTLVNegElectron,fPreviousEventTLVPosElectron);\r
1443     fPreviousEventTLVNegElectron = vCurrentTLVeNeg;\r
1444     fPreviousEventTLVPosElectron = vCurrentTLVePos;\r
1445 \r
1446   }\r
1447 \r
1448   /*\r
1449   //Photons P\r
1450   Double_t vtx[3];\r
1451   vtx[0]=0;vtx[1]=0;vtx[2]=0;\r
1452   for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){\r
1453 \r
1454   //      if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;\r
1455 \r
1456 \r
1457 \r
1458   Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();\r
1459   //      cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;\r
1460 \r
1461   //      cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;\r
1462 \r
1463   if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )\r
1464 \r
1465   fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());\r
1466 \r
1467 \r
1468   }\r
1469 \r
1470 \r
1471   */\r
1472 \r
1473 \r
1474 }\r
1475 \r
1476 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){\r
1477   //see header file for documentation\r
1478   for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){\r
1479     for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){\r
1480       fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));\r
1481     }\r
1482   }\r
1483 }\r
1484 void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, vector <TLorentzVector> eNeg, vector <TLorentzVector> ePos){\r
1485   //see header file for documentation\r
1486   for( UInt_t n=0; n < eNeg.size(); n++){\r
1487 \r
1488     TLorentzVector en = eNeg.at(n);\r
1489     for (UInt_t p=0; p < ePos.size(); p++){\r
1490       TLorentzVector ep = ePos.at(p);\r
1491       TLorentzVector np = ep + en;\r
1492       fHistograms->FillHistogram(histoName.Data(),np.M());\r
1493     }\r
1494   }\r
1495 \r
1496 }\r
1497 \r
1498 void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,vector <AliKFParticle> fKFGammas,\r
1499                                                               vector <TLorentzVector> tlVeNeg,vector<TLorentzVector> tlVePos)\r
1500 {\r
1501   //see header file for documentation\r
1502 \r
1503   for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++ ){\r
1504 \r
1505     for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){\r
1506 \r
1507       TLorentzVector xy = tlVePos[iPos] + tlVeNeg[iNeg];\r
1508 \r
1509       for (UInt_t iGam=0; iGam < fKFGammas.size(); iGam++){\r
1510 \r
1511         AliKFParticle * gammaCandidate = &fKFGammas[iGam];\r
1512         TLorentzVector g;\r
1513 \r
1514         g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);\r
1515         TLorentzVector xyg = xy + g;\r
1516         fHistograms->FillHistogram(histoMass.Data(),xyg.M());\r
1517         fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));\r
1518       }\r
1519     }\r
1520   }\r
1521 \r
1522 }\r
1523 void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, vector <TLorentzVector> e)\r
1524 {\r
1525   // see header file for documentation\r
1526   for(UInt_t i=0; i < e.size(); i++)\r
1527     {\r
1528       for (UInt_t j=i+1; j < e.size(); j++)\r
1529         {\r
1530           TLorentzVector ee = e[i] + e[j];\r
1531 \r
1532           fHistograms->FillHistogram(hBg.Data(),ee.M());\r
1533         }\r
1534     }\r
1535 }\r
1536 \r
1537 \r
1538 void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(vector <AliESDtrack*> negativeElectrons,\r
1539                                                         vector <AliESDtrack*> positiveElectrons, vector <AliKFParticle> gammas){\r
1540   // see header file for documentation\r
1541 \r
1542   UInt_t  sizeN = negativeElectrons.size();\r
1543   UInt_t  sizeP = positiveElectrons.size();\r
1544   UInt_t  sizeG = gammas.size();\r
1545 \r
1546 \r
1547 \r
1548   vector <Bool_t> xNegBand(sizeN);\r
1549   vector <Bool_t> xPosBand(sizeP);\r
1550   vector <Bool_t> gammaBand(sizeG);\r
1551 \r
1552 \r
1553   for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;\r
1554   for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;\r
1555   for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;\r
1556 \r
1557 \r
1558   for(UInt_t iPos=0; iPos < sizeP; iPos++){\r
1559         \r
1560     Double_t aP[3]; positiveElectrons[iPos]->GetConstrainedPxPyPz(aP); \r
1561 \r
1562     TVector3 ePosV(aP[0],aP[1],aP[2]);\r
1563 \r
1564     for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){\r
1565         \r
1566       Double_t aN[3]; negativeElectrons[iNeg]->GetConstrainedPxPyPz(aN); \r
1567       TVector3 eNegV(aN[0],aN[1],aN[2]);\r
1568 \r
1569       if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma\r
1570         xPosBand[iPos]=kFALSE;\r
1571         xNegBand[iNeg]=kFALSE;\r
1572       }\r
1573 \r
1574       for(UInt_t iGam=0; iGam < sizeG; iGam++){\r
1575         AliKFParticle* gammaCandidate = &gammas[iGam];\r
1576         TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());\r
1577         if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)\r
1578           gammaBand[iGam]=kFALSE;\r
1579       }\r
1580     }\r
1581   }\r
1582 \r
1583 \r
1584 \r
1585 \r
1586   for(UInt_t iPos=0; iPos < sizeP; iPos++){\r
1587     if(xPosBand[iPos]){\r
1588       fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);\r
1589     }\r
1590   }\r
1591   for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){\r
1592     if(xNegBand[iNeg]){\r
1593       fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);\r
1594     }\r
1595   }\r
1596   for(UInt_t iGam=0; iGam < sizeG; iGam++){\r
1597     if(gammaBand[iGam]){\r
1598       fKFReconstructedGammasCut.push_back(gammas[iGam]);\r
1599     }\r
1600   }\r
1601 }\r
1602 \r
1603 \r
1604 void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)\r
1605 {\r
1606   // see header file for documentation\r
1607   pid = -1;\r
1608   weight = -1;\r
1609 \r
1610   double wpart[5];\r
1611   double wpartbayes[5];\r
1612 \r
1613   //get probability of the diffenrent particle types\r
1614   track->GetESDpid(wpart);\r
1615 \r
1616   // Tentative particle type "concentrations"\r
1617   double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};\r
1618 \r
1619   //Bayes' formula\r
1620   double rcc = 0.;\r
1621   for (int i = 0; i < 5; i++)\r
1622     {\r
1623       rcc+=(c[i] * wpart[i]);\r
1624     }\r
1625 \r
1626 \r
1627 \r
1628   for (int i=0; i<5; i++) {\r
1629     if( rcc!=0){\r
1630       wpartbayes[i] = c[i] * wpart[i] / rcc;\r
1631     }\r
1632   }\r
1633 \r
1634 \r
1635 \r
1636   Float_t max=0.;\r
1637   int ipid=-1;\r
1638   //find most probable particle in ESD pid\r
1639   //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton\r
1640   for (int i = 0; i < 5; i++)\r
1641     {\r
1642       if (wpartbayes[i] > max)\r
1643         {\r
1644           ipid = i;\r
1645           max = wpartbayes[i];\r
1646         }\r
1647     }\r
1648 \r
1649   pid = ipid;\r
1650   weight = max;\r
1651 }\r
1652 double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)\r
1653 {\r
1654   // Calculates the number of sigma to the vertex.\r
1655 \r
1656   Float_t b[2];\r
1657   Float_t bRes[2];\r
1658   Float_t bCov[3];\r
1659   t->GetImpactParameters(b,bCov);\r
1660   if (bCov[0]<=0 || bCov[2]<=0) {\r
1661     AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1662     bCov[0]=0; bCov[2]=0;\r
1663   }\r
1664   bRes[0] = TMath::Sqrt(bCov[0]);\r
1665   bRes[1] = TMath::Sqrt(bCov[2]);\r
1666 \r
1667   // -----------------------------------\r
1668   // How to get to a n-sigma cut?\r
1669   //\r
1670   // The accumulated statistics from 0 to d is\r
1671   //\r
1672   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)\r
1673   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)\r
1674   //\r
1675   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)\r
1676   // Can this be expressed in a different way?\r
1677 \r
1678   if (bRes[0] == 0 || bRes[1] ==0)\r
1679     return -1;\r
1680 \r
1681   double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));\r
1682 \r
1683   // stupid rounding problem screws up everything:\r
1684   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(\r
1685   if (TMath::Exp(-d * d / 2) < 1e-10)\r
1686     return 1000;\r
1687 \r
1688 \r
1689   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);\r
1690   return d;\r
1691 }\r
1692 vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){\r
1693 \r
1694   vector <TLorentzVector> tlVtrack(0);\r
1695 \r
1696   for(UInt_t itrack=0; itrack < esdTrack.size(); itrack++){\r
1697     double P[3]; esdTrack[itrack]->GetConstrainedPxPyPz(P);\r
1698     TLorentzVector currentTrack;\r
1699     currentTrack.SetXYZM(P[0],P[1],P[2],fElectronMass);\r
1700     tlVtrack.push_back(currentTrack);\r
1701   }\r
1702 \r
1703   return tlVtrack;\r
1704 }\r
1705 \r