sync
authordlohner <dlohner@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 9 Aug 2012 11:45:12 +0000 (11:45 +0000)
committerdlohner <dlohner@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 9 Aug 2012 11:45:12 +0000 (11:45 +0000)
PWGGA/GammaConv/AliConversionAODBGHandlerRP.cxx [new file with mode: 0644]
PWGGA/GammaConv/AliConversionAODBGHandlerRP.h [new file with mode: 0644]
PWGGA/GammaConv/AliV0ReaderV1.cxx [new file with mode: 0644]
PWGGA/GammaConv/AliV0ReaderV1.h [new file with mode: 0644]

diff --git a/PWGGA/GammaConv/AliConversionAODBGHandlerRP.cxx b/PWGGA/GammaConv/AliConversionAODBGHandlerRP.cxx
new file mode 100644 (file)
index 0000000..c04270a
--- /dev/null
@@ -0,0 +1,270 @@
+#if !defined( __CINT__) || defined(__MAKECINT__)
+
+#include <exception>
+#include <iostream>
+#include "AliLog.h"
+#include "AliConversionAODBGHandlerRP.h"
+using namespace std;
+#endif
+
+
+// Author Daniel Lohner (Daniel.Lohner@cern.ch)
+
+
+ClassImp(AliConversionAODBGHandlerRP);
+
+//________________________________________________________________________
+AliConversionAODBGHandlerRP::AliConversionAODBGHandlerRP(Bool_t IsHeavyIon,Bool_t UseChargedTrackMult,Int_t NEvents) : TObject(),
+    fIsHeavyIon(IsHeavyIon),
+    fUseChargedTrackMult(UseChargedTrackMult),
+    fNEvents(NEvents),
+    fBGEventCounter(NULL),
+    fNBGEvents(NULL),
+    fNBinsZ(7),
+    fNBinsMultiplicity(5+Int_t(fUseChargedTrackMult)),
+    fBinLimitsArrayZ(NULL),
+    fBinLimitsArrayMultiplicity(NULL),
+    fBGPool(fNBinsZ,AliGammaConversionMultiplicityVector(fNBinsMultiplicity,AliGammaConversionBGEventVector(fNEvents)))
+{
+    // Vertex Z Binning
+
+    fBinLimitsArrayZ = new Double_t[fNBinsZ+1];
+    fBinLimitsArrayZ[0] = -50.00;
+    fBinLimitsArrayZ[1] = -3.375;
+    fBinLimitsArrayZ[2] = -1.605;
+    fBinLimitsArrayZ[3] = -0.225;
+    fBinLimitsArrayZ[4] = 1.065;
+    fBinLimitsArrayZ[5] = 2.445;
+    fBinLimitsArrayZ[6] = 4.245;
+    fBinLimitsArrayZ[7] = 50.00;
+
+    // MultiplicityBins
+
+    fBinLimitsArrayMultiplicity= new Double_t[fNBinsMultiplicity+1];
+
+    if(fUseChargedTrackMult){
+        // Use Charged Particle Multiplicity
+
+       fBinLimitsArrayMultiplicity[0] = 0;
+       fBinLimitsArrayMultiplicity[1] = 8.5;
+       fBinLimitsArrayMultiplicity[2] = 16.5;
+       fBinLimitsArrayMultiplicity[3] = 27.5;
+       fBinLimitsArrayMultiplicity[4] = 41.5;
+       fBinLimitsArrayMultiplicity[5] = 200.;
+
+       if(fIsHeavyIon){
+           fBinLimitsArrayMultiplicity[0] = 0;
+           fBinLimitsArrayMultiplicity[1] = 200.;
+           fBinLimitsArrayMultiplicity[2] = 500.;
+           fBinLimitsArrayMultiplicity[3] = 1000.;
+           fBinLimitsArrayMultiplicity[4] = 1500.;
+           fBinLimitsArrayMultiplicity[5] = 5000.;
+       }
+    }
+    else{
+       // Use V0 Multiplicity
+
+       fBinLimitsArrayMultiplicity[0] = 2;
+       fBinLimitsArrayMultiplicity[1] = 3;
+       fBinLimitsArrayMultiplicity[2] = 4;
+       fBinLimitsArrayMultiplicity[3] = 5;
+       fBinLimitsArrayMultiplicity[4] = 9999;
+
+       if(fIsHeavyIon){
+           fBinLimitsArrayMultiplicity[0] = 2;
+           fBinLimitsArrayMultiplicity[1] = 10;
+           fBinLimitsArrayMultiplicity[2] = 30;
+           fBinLimitsArrayMultiplicity[3] = 50;
+           fBinLimitsArrayMultiplicity[4] = 9999;
+       }
+    }
+    Initialize();
+}
+
+//________________________________________________________________________
+AliConversionAODBGHandlerRP::~AliConversionAODBGHandlerRP()
+{
+    if(fBinLimitsArrayZ){
+       delete[] fBinLimitsArrayZ;
+       fBinLimitsArrayZ=0x0;}
+    if(fBinLimitsArrayMultiplicity){
+       delete[] fBinLimitsArrayMultiplicity;
+       fBinLimitsArrayMultiplicity=0x0;}
+
+    if(fBGEventCounter){
+
+       for(Int_t z=0;z<fNBinsZ;z++){
+           delete[] fBGEventCounter[z];
+       }
+       delete[] fBGEventCounter;
+       fBGEventCounter = NULL;
+    }
+
+    // Delete pool
+
+    for(Int_t z=0;z<fNBinsZ;z++){
+       for(Int_t m=0;m<fNBinsMultiplicity;m++){
+           for(Int_t eventCounter=0;eventCounter<fNBGEvents[z][m]&&eventCounter<fNEvents;eventCounter++){
+
+               for(UInt_t d=0;d<fBGPool[z][m][eventCounter].size();d++){
+                   delete (AliAODConversionPhoton*)(fBGPool[z][m][eventCounter][d]);
+               }
+           }
+       }
+    }
+
+    if(fNBGEvents){
+       for(Int_t z=0;z<fNBinsZ;z++){
+           delete[] fNBGEvents[z];
+       }
+       delete[] fNBGEvents;
+       fNBGEvents = NULL;
+    }
+
+}
+
+//________________________________________________________________________
+void AliConversionAODBGHandlerRP::Initialize(){
+
+    // Counter
+
+    if(fBGEventCounter == NULL){
+       fBGEventCounter= new Int_t*[fNBinsZ];
+    }
+    for(Int_t z=0;z<fNBinsZ;z++){
+       fBGEventCounter[z]=new Int_t[fNBinsMultiplicity];
+    }
+
+    for(Int_t z=0;z<fNBinsZ;z++){
+       for(Int_t m=0;m<fNBinsMultiplicity;m++){
+           fBGEventCounter[z][m]=0;
+       }
+    }
+
+    if(fNBGEvents == NULL){
+       fNBGEvents= new Int_t*[fNBinsZ];
+    }
+    for(Int_t z=0;z<fNBinsZ;z++){
+       fNBGEvents[z]=new Int_t[fNBinsMultiplicity];
+    }
+    for(Int_t z=0;z<fNBinsZ;z++){
+       for(Int_t m=0;m<fNBinsMultiplicity;m++){
+           fNBGEvents[z][m]=0;
+       }
+    }
+}
+
+//-------------------------------------------------------------
+
+Int_t AliConversionAODBGHandlerRP::GetZBinIndex(Double_t zvalue) const{
+
+    if(zvalue<=fBinLimitsArrayZ[0]){
+       return -1;
+    }
+
+    if(fNBinsZ<2){return 0;}
+
+    for(Int_t i=0; i<fNBinsZ ;i++){
+       if(zvalue >= fBinLimitsArrayZ[i] && zvalue <= fBinLimitsArrayZ[i+1]){
+           return i;
+       }
+    }
+    return -1;
+}
+//-------------------------------------------------------------
+Int_t AliConversionAODBGHandlerRP::GetMultiplicityBinIndex(Int_t multiplicity) const{
+  if(fNBinsMultiplicity<2){
+    return 0;
+  }
+
+  for(Int_t i=0; i<fNBinsMultiplicity ;i++){
+    if(multiplicity >= fBinLimitsArrayMultiplicity[i] && multiplicity < fBinLimitsArrayMultiplicity[i+1]){
+      return i;
+    }
+  }
+  return -1;
+}
+
+//-------------------------------------------------------------
+Bool_t AliConversionAODBGHandlerRP::FindBins(TObjArray * const eventGammas,AliVEvent *fInputEvent,Int_t &zbin,Int_t &mbin){
+    Double_t vertexz=fInputEvent->GetPrimaryVertex()->GetZ();
+    zbin = GetZBinIndex(vertexz);
+
+    Int_t multiplicity=0;
+    if(fUseChargedTrackMult){
+       multiplicity=fInputEvent->GetNumberOfTracks();
+    }
+    else{
+       multiplicity=eventGammas->GetEntries();
+    }
+    mbin = GetMultiplicityBinIndex(multiplicity);
+
+    if(zbin<fNBinsZ&&mbin<fNBinsMultiplicity){
+       if(zbin>=0&&mbin>=0){
+           return kTRUE;
+       }
+    }
+    //cout<<Form("Requested BG pool does not exist:  z %i m %i",zbin,mbin)<<endl;
+    return kFALSE;
+}
+
+
+//-------------------------------------------------------------
+void AliConversionAODBGHandlerRP::AddEvent(TObjArray * const eventGammas,AliVEvent *fInputEvent){
+
+    if(eventGammas->GetEntriesFast()==0)return;
+
+    Int_t z;
+    Int_t m;
+    if(FindBins(eventGammas,fInputEvent,z,m)){
+  
+   
+       // If Event Stack is full, replace the first entry (First in first out)
+       if(fBGEventCounter[z][m] >= fNEvents){
+           fBGEventCounter[z][m]=0;
+       }
+
+       // Update number of Events stored
+       if(fNBGEvents[z][m] < fNEvents){
+           fNBGEvents[z][m]++;
+       }
+
+       Int_t eventCounter=fBGEventCounter[z][m];
+
+       //clear the vector for old gammas
+       for(UInt_t d=0;d<fBGPool[z][m][eventCounter].size();d++){
+           delete (AliAODConversionPhoton*)(fBGPool[z][m][eventCounter][d]);
+       }
+
+       fBGPool[z][m][eventCounter].clear();
+
+       // add the gammas to the vector
+
+       for(Int_t i=0; i< eventGammas->GetEntriesFast();i++){
+           fBGPool[z][m][eventCounter].push_back(new AliAODConversionPhoton(*(AliAODConversionPhoton*)(eventGammas->At(i))));
+       }
+
+       fBGEventCounter[z][m]++;
+    }
+}
+
+//-------------------------------------------------------------
+AliGammaConversionPhotonVector* AliConversionAODBGHandlerRP::GetBGGoodGammas(TObjArray * const eventGammas,AliVEvent *fInputEvent,Int_t event){
+    Int_t zbin;
+    Int_t mbin;
+    if(FindBins(eventGammas,fInputEvent,zbin,mbin)){
+       return &(fBGPool[zbin][mbin][event]);
+    }
+    return NULL;
+}
+
+//-------------------------------------------------------------
+Int_t AliConversionAODBGHandlerRP::GetNBGEvents(TObjArray * const eventGammas,AliVEvent *fInputEvent){
+    Int_t zbin;
+    Int_t mbin;
+    if(FindBins(eventGammas,fInputEvent,zbin,mbin)){
+       return fNBGEvents[zbin][mbin];
+    }
+    return 0;
+}
+
diff --git a/PWGGA/GammaConv/AliConversionAODBGHandlerRP.h b/PWGGA/GammaConv/AliConversionAODBGHandlerRP.h
new file mode 100644 (file)
index 0000000..144eb98
--- /dev/null
@@ -0,0 +1,52 @@
+#ifndef __ALICONVERSIONAODBGHANDLERRP_H__
+#define __ALICONVERSIONAODBGHANDLERRP_H__
+
+#include "AliLog.h"
+#include "TObject.h"
+#include "AliAODConversionPhoton.h"
+#include "TObjArray.h"
+using namespace std;
+
+typedef vector<AliAODConversionPhoton*> AliGammaConversionPhotonVector;   // Vector containing photons
+typedef vector<AliGammaConversionPhotonVector> AliGammaConversionBGEventVector;       // Event contains vector of gammas (AliConversionPhotons)
+typedef vector<AliGammaConversionBGEventVector> AliGammaConversionMultiplicityVector;  // Multiplicity classes containing event vectors
+typedef vector<AliGammaConversionMultiplicityVector> AliGammaConversionBGVector;       // z vertex position ...
+
+
+class AliConversionAODBGHandlerRP: public TObject{
+
+public:
+
+    AliConversionAODBGHandlerRP(Bool_t IsHeavyIon=kFALSE,Bool_t UseChargedTrackMult=kTRUE,Int_t NEvents=10);
+    
+    virtual ~AliConversionAODBGHandlerRP();
+
+    Int_t GetZBinIndex(Double_t z) const;
+    Int_t GetMultiplicityBinIndex(Int_t mult) const;
+    void Initialize();
+    Bool_t FindBins(TObjArray * const eventGammas,AliVEvent *fInputEvent,Int_t &zbin,Int_t &mbin);
+
+    AliGammaConversionPhotonVector* GetBGGoodGammas(TObjArray * const eventGammas,AliVEvent *fInputEvent,Int_t event);
+    void AddEvent(TObjArray * const eventGammas,AliVEvent *fInputEvent);
+    Int_t GetNBGEvents()const {return fNEvents;} // Size of the Pool (20)
+    Int_t GetNBGEvents(TObjArray * const eventGammas,AliVEvent *fInputEvent);
+    Int_t GetNZBins()const{return fNBinsZ;};
+    Int_t GetNMultiplicityBins()const{return fNBinsMultiplicity;};
+
+private:
+    Bool_t fIsHeavyIon;
+    Bool_t fUseChargedTrackMult;
+    Int_t fNEvents;
+    Int_t **fBGEventCounter; //! bg counter
+    Int_t **fNBGEvents;
+    Int_t fNBinsZ; //n z bins
+    Int_t fNBinsMultiplicity; //n bins multiplicity
+    Double_t *fBinLimitsArrayZ;//! bin limits z array
+    Double_t *fBinLimitsArrayMultiplicity;//! bin limit multiplicity array
+    AliGammaConversionBGVector fBGPool; //background events
+
+ ClassDef(AliConversionAODBGHandlerRP,0);
+
+};
+#endif
diff --git a/PWGGA/GammaConv/AliV0ReaderV1.cxx b/PWGGA/GammaConv/AliV0ReaderV1.cxx
new file mode 100644 (file)
index 0000000..456ced0
--- /dev/null
@@ -0,0 +1,620 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                                                                                                                                                                                                                                             *
+ * Authors: Svein Lindal, Daniel Lohner                                                                                                *
+ * Version 1.0                                                                                                                                                                                                                                         *
+ *                                                                                                                                                                                                                                                                                             *
+ * Permission to use, copy, modify and distribute this software and its         *
+ * documentation strictly for non-commercial purposes is hereby granted         *
+ * without fee, provided that the above copyright notice appears in all         *
+ * copies and that both the copyright notice and this permission notice         *
+ * appear in the supporting documentation. The authors make no claims           *
+ * about the suitability of this software for any purpose. It is                                       *
+ * provided "as is" without express or implied warranty.                                                                       *
+ **************************************************************************/
+
+////////////////////////////////////////////////
+//--------------------------------------------- 
+// Class reconstructing conversion photons from V0s
+//---------------------------------------------
+////////////////////////////////////////////////
+
+#include "AliV0ReaderV1.h"
+#include "AliKFParticle.h"
+#include "AliAODv0.h"
+#include "AliESDv0.h"
+#include "AliAODEvent.h"
+#include "AliESDEvent.h"
+#include "AliKFParticle.h"
+#include "AliKFConversionPhoton.h"
+#include "AliAODConversionPhoton.h"
+#include "AliConversionPhotonBase.h"
+#include "TVector.h"
+#include "AliKFVertex.h"
+#include "AliAODTrack.h"
+#include "AliESDtrack.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+#include "AliAODHandler.h"
+#include "AliPIDResponse.h"
+#include "TChain.h"
+#include "AliStack.h"
+
+class iostream;
+
+
+using namespace std;
+
+ClassImp(AliV0ReaderV1)
+
+//________________________________________________________________________
+AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
+    fConversionCuts(NULL),
+    fConversionGammas(NULL),
+    fUseImprovedVertex(kTRUE),
+    fUseOwnXYZCalculation(kTRUE),
+    fUseConstructGamma(kFALSE),
+    kUseAODConversionPhoton(kTRUE),
+    fCreateAOD(kFALSE),
+    fDeltaAODBranchName("GammaConv"),
+    fDeltaAODFilename("AliAODGammaConversion.root"),
+    fEventIsSelected(kFALSE)
+{
+    // Default constructor
+
+    DefineInput(0, TChain::Class());
+}
+
+//________________________________________________________________________
+AliV0ReaderV1::~AliV0ReaderV1()
+{
+    // default deconstructor
+
+    if(fConversionGammas){
+       fConversionGammas->Delete();// Clear Objects
+       delete fConversionGammas;
+       fConversionGammas=0x0;
+    }
+}
+//________________________________________________________________________
+void AliV0ReaderV1::Init()
+{
+    // Initialize function to be called once before analysis
+
+    if(fConversionCuts==NULL){
+       if(fConversionCuts==NULL)AliError("No Cut Selection initialized");
+    }
+
+    if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
+
+    if(fConversionGammas != NULL){
+       delete fConversionGammas;
+       fConversionGammas=NULL;
+    }
+
+    if(fConversionGammas == NULL){
+       if(kUseAODConversionPhoton){
+           fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
+       else{
+           fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
+    }
+    fConversionGammas->Delete();//Reset the TClonesArray
+   
+    // Create AODs
+
+    if(fCreateAOD){
+       if(fConversionCuts){
+           fDeltaAODBranchName.Append("_");
+           fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
+           fDeltaAODBranchName.Append("_gamma");
+       }
+       fConversionGammas->SetName(fDeltaAODBranchName.Data());
+
+       AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
+       AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
+    }
+}
+
+//________________________________________________________________________
+void AliV0ReaderV1::UserCreateOutputObjects()
+{
+    // Create User Output Objects
+}
+
+//________________________________________________________________________
+void AliV0ReaderV1::UserExec(Option_t *){
+
+    // User Exec
+    fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
+}
+
+//________________________________________________________________________
+Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
+{
+    //Reset the TClonesArray
+    fConversionGammas->Delete();
+
+    fInputEvent=inputEvent;
+    fMCEvent=mcEvent;
+
+    if(!fInputEvent){
+       AliError("No Input event");
+       return kFALSE;
+    }
+
+    if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
+   
+    // Event Cuts
+    if(!fConversionCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
+
+    // Set Magnetic Field
+    AliKFParticle::SetField(fInputEvent->GetMagneticField());
+
+    if(fInputEvent->IsA()==AliESDEvent::Class()){
+       ProcessESDV0s();
+    }
+    if(fInputEvent->IsA()==AliAODEvent::Class()){
+       GetAODConversionGammas();
+    }
+
+    // AOD Output
+    FillAODOutput();
+
+    return kTRUE;
+}
+///________________________________________________________________________
+void AliV0ReaderV1::FillAODOutput()
+{
+    // Fill AOD Output with reconstructed Photons
+
+    if(fInputEvent->IsA()==AliESDEvent::Class()){
+       ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
+       if(fCreateAOD) {
+           AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+           if (aodhandler && aodhandler->GetFillAOD()) {
+             AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
+             PostData(0, fConversionGammas);
+             
+           }
+       }
+    }
+}
+
+///________________________________________________________________________
+const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
+
+    // Get External Track Parameter with given charge
+
+    if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
+
+    // Check for sign flip
+    if(fCurrentV0){
+       if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
+       if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
+       if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
+           tracklabel=fCurrentV0->GetPindex();
+           return fCurrentV0->GetParamP();}
+       if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
+           tracklabel=fCurrentV0->GetNindex();
+           return fCurrentV0->GetParamN();}
+    }
+    return 0x0;
+}
+
+///________________________________________________________________________
+Bool_t AliV0ReaderV1::ProcessESDV0s()
+{
+    // Process ESD V0s for conversion photon reconstruction
+
+    AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
+
+    AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
+   
+    if(fESDEvent){
+
+       for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
+           AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
+           if(!fCurrentV0){
+               printf("Requested V0 does not exist");
+               continue;}
+
+           fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
+
+           if(fCurrentMotherKFCandidate){
+
+               // Add Gamma to the TClonesArray
+
+               if(kUseAODConversionPhoton){
+                   new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
+               }
+               else{
+                   new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
+               }
+
+               delete fCurrentMotherKFCandidate;
+               fCurrentMotherKFCandidate=NULL;
+           }
+       }
+    }
+    return kTRUE;
+}
+
+///________________________________________________________________________
+AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
+{
+       // Reconstruct conversion photon from ESD v0
+    fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonIn);
+
+    //checks if on the fly mode is set
+    if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
+               fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kOnFly);
+               return 0x0;
+    }
+
+    // TrackLabels
+    Int_t currentTrackLabels[2]={-1,-1};
+
+    // Get Daughter KF Particles
+
+    const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
+    const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
+
+    if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
+
+    // Apply some Cuts before Reconstruction
+
+    AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
+    AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
+
+    if(!negTrack || !posTrack) {
+               fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kNoTracks);
+               return 0x0;
+    }
+
+    // Track Cuts
+    if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
+               fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kTrackCuts);
+               return 0x0;
+    }
+
+    // PID Cuts
+       if(!fConversionCuts->dEdxCuts(negTrack) || !fConversionCuts->dEdxCuts(posTrack)) {
+               fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
+               return 0x0;
+       }
+
+    // Reconstruct Photon
+
+    AliKFConversionPhoton *fCurrentMotherKF=NULL;
+
+    AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
+    AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
+
+    // Reconstruct Gamma
+
+    if(fUseConstructGamma){
+
+       fCurrentMotherKF = new AliKFConversionPhoton();
+       fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
+    }else{
+       fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
+       fCurrentMotherKF->SetMassConstraint(0,0);
+    }
+
+    // Set Track Labels
+
+    fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
+
+    // Set V0 index
+
+    fCurrentMotherKF->SetV0Index(currentV0Index);
+
+    //Set MC Label
+
+    if(fMCEvent){
+
+       AliStack *fMCStack= fMCEvent->Stack();
+
+       Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
+       Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
+
+       TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
+       TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
+
+       if(fPositiveMCParticle&&fNegativeMCParticle){
+           fCurrentMotherKF->SetMCLabelPositive(labelp);
+           fCurrentMotherKF->SetMCLabelNegative(labeln);
+       }
+    }
+
+    // Update Vertex (moved for same eta compared to old)
+    if(fUseImprovedVertex == kTRUE){
+       AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
+       primaryVertexImproved+=*fCurrentMotherKF;
+       fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
+    }
+
+    // SetPsiPair
+
+    Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
+    fCurrentMotherKF->SetPsiPair(PsiPair);
+
+    
+    // Recalculate ConversionPoint
+    Double_t dca[2]={0,0};
+    if(fUseOwnXYZCalculation){
+       Double_t convpos[3]={0,0,0};
+       if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
+           fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kConvPointFail);
+           delete fCurrentMotherKF;
+           fCurrentMotherKF=NULL;
+           return 0x0;
+       }
+
+       fCurrentMotherKF->SetConversionPoint(convpos);
+    }
+
+    if(fCurrentMotherKF->GetNDF() > 0.)
+       fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF());   //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
+    
+    // Apply Photon Cuts
+
+    if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
+       fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonCuts);
+       delete fCurrentMotherKF;
+       fCurrentMotherKF=NULL;
+       return 0x0;
+    }
+
+    // Set Dilepton Mass (moved down for same eta compared to old)
+    fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
+
+
+    fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonOut);
+    return fCurrentMotherKF;
+}
+
+///________________________________________________________________________
+Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
+    //
+    // Angle between daughter momentum plane and plane
+    //
+
+    AliExternalTrackParam nt(*negativeparam);
+    AliExternalTrackParam pt(*positiveparam);
+
+   Float_t magField = fInputEvent->GetMagneticField();
+
+   Double_t xyz[3] = {0.,0.,0.};
+   v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
+     
+   Double_t mn[3] = {0,0,0};
+   Double_t mp[3] = {0,0,0};
+  
+   v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
+   v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter; 
+
+   Double_t deltat = 1.;
+   deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) -  TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
+   Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated 
+   Double_t momPosProp[3] = {0,0,0};
+   Double_t momNegProp[3] = {0,0,0};
+    
+   Double_t psiPair = 4.;
+   if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
+   
+   if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
+  
+   pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
+   nt.GetPxPyPz(momNegProp);
+
+   Double_t pEle =
+       TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
+
+   Double_t pPos =
+       TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
+
+   Double_t scalarproduct =
+       momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
+
+   Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
+
+   psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
+
+   return psiPair;
+}
+
+///________________________________________________________________________
+Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
+
+    // Get Center of the helix track parametrization
+
+    Int_t charge=track->Charge();
+    Double_t b=fInputEvent->GetMagneticField();
+
+    Double_t   helix[6];
+    track->GetHelixParameters(helix,b);
+
+    Double_t xpos =    helix[5];
+    Double_t ypos =    helix[0];
+    Double_t radius = TMath::Abs(1./helix[4]);
+    Double_t phi = helix[2];
+
+    if(phi < 0){
+       phi = phi + 2*TMath::Pi();
+    }
+
+    phi -= TMath::Pi()/2.;
+    Double_t xpoint =  radius * TMath::Cos(phi);
+    Double_t ypoint =  radius * TMath::Sin(phi);
+
+    if(b<0){
+       if(charge > 0){
+           xpoint = - xpoint;
+           ypoint = - ypoint;
+       }
+
+       if(charge < 0){
+           xpoint =    xpoint;
+           ypoint =    ypoint;
+       }
+    }
+    if(b>0){
+       if(charge > 0){
+           xpoint =    xpoint;
+           ypoint =    ypoint;
+       }
+
+       if(charge < 0){
+           xpoint = - xpoint;
+           ypoint = - ypoint;
+       }
+    }
+    center[0] =        xpos + xpoint;
+    center[1] =        ypos + ypoint;
+
+    return 1;
+}
+///________________________________________________________________________
+Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
+
+    // Recalculate Conversion Point
+
+    if(!pparam||!nparam)return kFALSE;
+    Double_t helixcenterpos[2];
+    GetHelixCenter(pparam,helixcenterpos);
+
+    Double_t helixcenterneg[2];
+    GetHelixCenter(nparam,helixcenterneg);
+
+       Double_t helixpos[6];
+       pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
+       Double_t posradius = TMath::Abs(1./helixpos[4]);
+
+       Double_t helixneg[6];
+       nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
+       Double_t negradius = TMath::Abs(1./helixneg[4]);
+
+        // Calculate xy-position
+
+       Double_t xpos = helixcenterpos[0];
+       Double_t ypos = helixcenterpos[1];
+       Double_t xneg = helixcenterneg[0];
+       Double_t yneg = helixcenterneg[1];
+
+       convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
+       convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
+
+
+       // Calculate Track XY vertex position
+
+       Double_t deltaXPos = convpos[0] -       xpos;
+       Double_t deltaYPos = convpos[1] -       ypos;
+
+       Double_t deltaXNeg = convpos[0] -       xneg;
+       Double_t deltaYNeg = convpos[1] -       yneg;
+
+       Double_t alphaPos =     TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
+       Double_t alphaNeg =     TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
+
+       Double_t vertexXNeg =   xneg +  TMath::Abs(negradius)*TMath::Cos(alphaNeg);
+       Double_t vertexYNeg =   yneg +  TMath::Abs(negradius)*TMath::Sin(alphaNeg);
+
+       Double_t vertexXPos =   xpos +  TMath::Abs(posradius)*TMath::Cos(alphaPos);
+       Double_t vertexYPos =   ypos +  TMath::Abs(posradius)*TMath::Sin(alphaPos);
+
+       AliExternalTrackParam p(*pparam);
+        AliExternalTrackParam n(*nparam);
+
+       TVector2 vertexPos(vertexXPos,vertexYPos);
+       TVector2 vertexNeg(vertexXNeg,vertexYNeg);
+
+       // Convert to local coordinate system
+       TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
+       TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
+
+       // Propagate Track Params to Vertex
+
+       if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
+       if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
+
+       // Check whether propagation was sucessful
+
+       if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
+       if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
+
+       // Calculate z position
+
+       convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
+
+       // Calculate DCA
+       TVector2 vdca=vertexPos-vertexNeg;
+       dca[0]=vdca.Mod();
+       dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
+
+       return kTRUE;
+}
+//________________________________________________________________________
+Bool_t AliV0ReaderV1::GetAODConversionGammas(){
+
+    // Get reconstructed conversion photons from satellite AOD file
+
+    AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
+
+    if(fAODEvent){
+
+       if(fConversionGammas == NULL){
+           fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
+       }
+       fConversionGammas->Delete();//Reset the TClonesArray
+
+       //Get Gammas from satellite AOD gamma branch
+
+       AliAODConversionPhoton *gamma=0x0;
+    
+       TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
+       if(!fInputGammas){
+           FindDeltaAODBranchName();
+           fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
+       if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
+       // Apply Selection Cuts to Gammas and create local working copy
+       if(fInputGammas){
+           for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
+               gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
+               if(gamma){
+                   if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
+                       new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
+               }
+           }
+       }
+    }
+
+    if(fConversionGammas->GetEntries()){return kTRUE;}
+
+    return kFALSE;
+}
+
+//________________________________________________________________________
+void AliV0ReaderV1::FindDeltaAODBranchName(){
+
+    // Find delta AOD branchname containing reconstructed photons
+
+    TList *list=fInputEvent->GetList();
+    for(Int_t ii=0;ii<list->GetEntries();ii++){
+       TString name((list->At(ii))->GetName());
+       if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
+           fDeltaAODBranchName=name;
+           AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
+           return;
+       }
+    }
+}
+
+
+//________________________________________________________________________
+void AliV0ReaderV1::Terminate(Option_t *)
+{
+
+}
diff --git a/PWGGA/GammaConv/AliV0ReaderV1.h b/PWGGA/GammaConv/AliV0ReaderV1.h
new file mode 100644 (file)
index 0000000..c110cff
--- /dev/null
@@ -0,0 +1,103 @@
+#ifndef ALIV0READERV1_H
+#define ALIV0READERV1_H
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAODv0.h"
+#include "AliESDv0.h"
+#include "AliConversionCuts.h"
+#include "AliExternalTrackParam.h"
+
+class AliConversionPhotonBase;
+class TRandom3;
+class AliStack;
+class TList;
+class AliKFConversionPhoton;
+class TString;
+class TClonesArray;
+class TH1F;
+class TH2F;
+class AliAODConversionPhoton;
+
+using namespace std;
+
+class AliV0ReaderV1 : public AliAnalysisTaskSE {
+       
+ public: 
+       
+  AliV0ReaderV1(const char *name="V0ReaderV1");
+  virtual ~AliV0ReaderV1();                            //virtual destructor
+  void UserCreateOutputObjects();
+
+  virtual void UserExec(Option_t *option);
+  virtual void Terminate(Option_t *);
+  virtual void Init();
+
+  Bool_t ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent=NULL);
+  Bool_t IsEventSelected(){return fEventIsSelected;}
+
+  // Return Reconstructed Gammas
+  TClonesArray *GetReconstructedGammas(){return fConversionGammas;}
+  Int_t GetNReconstructedGammas(){if(fConversionGammas){return fConversionGammas->GetEntriesFast();}else{return 0;}}
+  AliConversionCuts *GetConversionCuts(){return fConversionCuts;}
+  TList *GetCutHistograms(){if(fConversionCuts){return fConversionCuts->GetCutHistograms();}return NULL;}
+  // Set Options
+
+  void SetConversionCuts(const TString cut);
+  void SetConversionCuts(AliConversionCuts *cuts){fConversionCuts=cuts;}
+  void SetUseOwnXYZCalculation(Bool_t flag){fUseOwnXYZCalculation=flag;}
+  void SetUseConstructGamma(Bool_t flag){fUseConstructGamma=flag;}
+  void SetUseAODConversionPhoton(Bool_t b){if(b){cout<<"Setting Outputformat to AliAODConversionPhoton "<<endl;}else{cout<<"Setting Outputformat to AliKFConversionPhoton "<<endl;};kUseAODConversionPhoton=b;}
+  void SetCreateAODs(Bool_t k=kTRUE){fCreateAOD=k;}
+  void SetDeltaAODFilename(TString s){fDeltaAODFilename=s;}
+  void SetDeltaAODBranchName(TString string) { fDeltaAODBranchName = string;AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));}
+
+protected:
+    // Reconstruct Gammas
+    Bool_t ProcessESDV0s();
+    AliKFConversionPhoton *ReconstructV0(AliESDv0* fCurrentV0,Int_t currentV0Index);
+    void FillAODOutput();
+    void FindDeltaAODBranchName();
+    Bool_t GetAODConversionGammas();
+   
+    // Getter Functions
+
+    const AliExternalTrackParam *GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge);
+    const AliExternalTrackParam *GetExternalTrackParamP(AliESDv0 *fCurrentV0,Int_t &tracklabel){return GetExternalTrackParam(fCurrentV0,tracklabel,1);};
+    const AliExternalTrackParam *GetExternalTrackParamN(AliESDv0 *fCurrentV0,Int_t &tracklabel){return GetExternalTrackParam(fCurrentV0,tracklabel,-1);};
+    AliKFParticle *GetPositiveKFParticle(AliAODv0 *fCurrentV0,Int_t fTrackLabel[2]);
+    AliKFParticle *GetNegativeKFParticle(AliAODv0 *fCurrentV0,Int_t fTrackLabel[2]);
+    AliKFParticle *GetPositiveKFParticle(AliESDv0 *fCurrentV0,Int_t fTrackLabel[2]);
+    AliKFParticle *GetNegativeKFParticle(AliESDv0 *fCurrentV0,Int_t fTrackLabel[2]);
+
+    Bool_t GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]);
+    Bool_t GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]);
+    Double_t GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const;
+
+    AliConversionCuts *fConversionCuts; // Pointer to the ConversionCut Selection
+    TClonesArray *fConversionGammas; // TClonesArray holding the reconstructed photons
+    Bool_t fUseImprovedVertex; // set flag to improve primary vertex estimation by adding photons
+    Bool_t fUseOwnXYZCalculation; //flag that determines if we use our own calculation of xyz (markus)
+    Bool_t fUseConstructGamma; //flag that determines if we use ConstructGamma method from AliKF
+    Bool_t kUseAODConversionPhoton; // set flag to use AOD instead of KF output format for photons
+    Bool_t fCreateAOD; // set flag for AOD creation
+    TString     fDeltaAODBranchName;// File where Gamma Conv AOD is located, if not in default AOD
+    TString fDeltaAODFilename; // set filename for delta/satellite aod
+    Bool_t fEventIsSelected;
+
+
+    ClassDef(AliV0ReaderV1,1)
+};
+
+inline void AliV0ReaderV1::SetConversionCuts(const TString cut){
+    if(fConversionCuts != NULL){
+       delete fConversionCuts;
+       fConversionCuts=NULL;
+    }
+    if(fConversionCuts == NULL){
+       fConversionCuts=new AliConversionCuts("V0ReaderCuts","V0ReaderCuts");
+       fConversionCuts->InitializeCutsFromCutString(cut.Data());
+    }
+}
+
+
+#endif