PWGGA/GammaConv/AliV0ReaderV1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Authors: Svein Lindal, Daniel Lohner                                   *
   5  * Version 1.0                                                            *
   6  *                                                                        *
   7  *                                                                        *
   8  * based on: on older version (see aliroot up to v5-04-42-AN)             *
   9  *           AliV0Reader.cxx                                              *
  10  *           Authors: Kathrin Koch, Kenneth Aamodt, Ana Marin             *
  11  *                                                                        *
  12  * Permission to use, copy, modify and distribute this software and its   *
  13  * documentation strictly for non-commercial purposes is hereby granted   *
  14  * without fee, provided that the above copyright notice appears in all   *
  15  * copies and that both the copyright notice and this permission notice   *
  16  * appear in the supporting documentation. The authors make no claims     *
  17  * about the suitability of this software for any purpose. It is          *
  18  * provided "as is" without express or implied warranty.                  *
  19  **************************************************************************/
  20
  21 ////////////////////////////////////////////////
  22 //---------------------------------------------
  23 // Class reconstructing conversion photons from V0s
  24 //---------------------------------------------
  25 ////////////////////////////////////////////////
  26
  27 // The AliAODConversionPhotons will return the Position (== ESDTrack->GetID())of the positive (negative) ESDTrack
  28 // in the Array of ESDTracks stored in the ESDEvent via GetTrackLabelPositive() (GetTrackLabelNagative()).
  29 // Since it is the Position in the Array it is always positive.
  30 // After the conversion ESD --> AOD each AOD track will give you the former Position in the ESDArray via the
  31 // Function AODTrack->GetID(). AODTracks are stored for different TrackParameter (e.g. HybridTracks, TPC only ...). No Standard
  32 // AODTracks (not copies of AliExternalTrackParam) will be stored with negative values of AODTrack->GetID() (Formula: (ESDTrack->GetID()+1)*-1)
  33 // This leads to the possibility of more than one AOD track with the same negative ID. For that reason all AOD tracks are additionally flaged.
  34 // In this analysis we should therefore only use AODTracks with positive values of AODTrack->GetID().
  35
  36 // If you want to find the AODTrack corresponding to the daugher track of a AliAODConversionPhoton you have to find the AODTrack with
  37 // AODTrack->GetID() == GetTrackLabelPositive() (GetTrackLabelNagative()).
  38
  39 #include <TGeoGlobalMagField.h>
  40
  41 #include "AliV0ReaderV1.h"
  42 #include "AliKFParticle.h"
  43 #include "AliAODv0.h"
  44 #include "AliESDv0.h"
  45 #include "AliAODEvent.h"
  46 #include "AliESDEvent.h"
  47 #include "AliPID.h"
  48 #include "AliMCEvent.h"
  49 #include "AliStack.h"
  50 #include "AliMCEventHandler.h"
  51 #include "AliESDpid.h"
  52 #include "AliESDtrackCuts.h"
  53 #include "TRandom3.h"
  54 #include "AliGenCocktailEventHeader.h"
  55 #include "TList.h"
  56 #include "AliKFConversionPhoton.h"
  57 #include "AliAODConversionPhoton.h"
  58 #include "AliConversionPhotonBase.h"
  59 #include "TVector.h"
  60 #include "AliKFVertex.h"
  61 #include "AliAODTrack.h"
  62 #include "AliESDtrack.h"
  63 #include "AliAnalysisManager.h"
  64 #include "AliInputEventHandler.h"
  65 #include "AliAODHandler.h"
  66 #include "AliPIDResponse.h"
  67 #include "TChain.h"
  68 #include "TFile.h"
  69 #include "TString.h"
  70 #include "TObjArray.h"
  71
  72 class iostream;
  73
  74
  75 using namespace std;
  76
  77 ClassImp(AliV0ReaderV1)
  78
  79 //________________________________________________________________________
  80 AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
  81    fConversionCuts(NULL),
  82    fConversionGammas(NULL),
  83    fUseImprovedVertex(kTRUE),
  84    fUseOwnXYZCalculation(kTRUE),
  85    fUseConstructGamma(kFALSE),
  86    kUseAODConversionPhoton(kTRUE),
  87    fCreateAOD(kFALSE),
  88    fDeltaAODBranchName("GammaConv"),
  89    fDeltaAODFilename("AliAODGammaConversion.root"),
  90    fRelabelAODs(kFALSE),
  91    fEventIsSelected(kFALSE),
  92    fNumberOfPrimaryTracks(0),
  93    fPeriodName("")
  94 {
  95    // Default constructor
  96
  97    DefineInput(0, TChain::Class());
  98 }
  99
 100 //________________________________________________________________________
 101 AliV0ReaderV1::~AliV0ReaderV1()
 102 {
 103    // default deconstructor
 104
 105    if(fConversionGammas){
 106       fConversionGammas->Delete();// Clear Objects
 107       delete fConversionGammas;
 108       fConversionGammas=0x0;
 109    }
 110 }
 111 /*
 112 //________________________________________________________________________
 113 AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
 114 fConversionCuts(NULL),
 115 fConversionGammas(NULL),
 116 fUseImprovedVertex(original.fUseImprovedVertex),
 117 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
 118 fUseConstructGamma(original.fUseConstructGamma),
 119 kUseAODConversionPhoton(original.kUseAODConversionPhoton),
 120 fCreateAOD(original.fCreateAOD),
 121 fDeltaAODBranchName(original.fDeltaAODBranchName),
 122 fDeltaAODFilename(original.fDeltaAODFilename),
 123 fEventIsSelected(original.fEventIsSelected)
 124 {
 125 // Default constructor
 126
 127 DefineInput(0, TChain::Class());
 128 }
 129
 130 //____________________________________________________________
 131 AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
 132 //
 133 // Assignment operator
 134 // Only copies pointers, object is not the owner of the references
 135 //
 136 if(this != &ref){
 137 AliAnalysisTaskSE::operator=(ref);
 138 fUseImprovedVertex=ref.fUseImprovedVertex;
 139 fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
 140 fUseConstructGamma=ref.fUseConstructGamma;
 141 kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
 142 fCreateAOD=ref.fCreateAOD;
 143 fDeltaAODBranchName=ref.fDeltaAODBranchName;
 144 fDeltaAODFilename=ref.fDeltaAODFilename;
 145 fEventIsSelected=ref.fEventIsSelected;
 146 }
 147 return *this;
 148 }
 149 */
 150 //________________________________________________________________________
 151 void AliV0ReaderV1::Init()
 152 {
 153    // Initialize function to be called once before analysis
 154    if(fConversionCuts==NULL){
 155       if(fConversionCuts==NULL)AliError("No Cut Selection initialized");
 156    }
 157    if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
 158
 159    if(fConversionGammas != NULL){
 160       delete fConversionGammas;
 161       fConversionGammas=NULL;
 162    }
 163
 164    if(fConversionGammas == NULL){
 165       if(kUseAODConversionPhoton){
 166          fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
 167       else{
 168          fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
 169    }
 170    fConversionGammas->Delete();//Reset the TClonesArray
 171 }
 172
 173 //________________________________________________________________________
 174 void AliV0ReaderV1::UserCreateOutputObjects()
 175 {
 176    // Create AODs
 177
 178    if(fCreateAOD){
 179       if(fConversionCuts){
 180          fDeltaAODBranchName.Append("_");
 181          fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
 182          fDeltaAODBranchName.Append("_gamma");
 183       }
 184       fConversionGammas->SetName(fDeltaAODBranchName.Data());
 185
 186       AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
 187       AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
 188    }
 189
 190 }
 191 //________________________________________________________________________
 192 Bool_t AliV0ReaderV1::Notify()
 193 {
 194    if (fPeriodName.CompareTo("") == 0){
 195       AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 196       if(man) {
 197          AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 198          if (inputHandler){
 199             TTree* tree = (TTree*) inputHandler->GetTree();
 200             TFile* file = (TFile*) tree->GetCurrentFile();
 201             TString fileName(file->GetName());
 202             TObjArray *arr = fileName.Tokenize("/");
 203             for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
 204                TObjString* testObjString = (TObjString*)arr->At(i);
 205                if (testObjString->GetString().Contains("LHC")){
 206                   fPeriodName = testObjString->GetString();
 207                   i = arr->GetEntriesFast();
 208                }
 209             }
 210          }
 211       }
 212       if(!fConversionCuts->GetDoEtaShift()) return kTRUE; // No Eta Shift requested, continue
 213       if(fConversionCuts->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
 214          fConversionCuts->GetCorrectEtaShiftFromPeriod(fPeriodName);
 215          fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 216          return kTRUE;
 217       }
 218       else{
 219          printf(" Gamma Conversion Reader %s :: Eta Shift Manually Set to %f \n\n",
 220                 (fConversionCuts->GetCutNumber()).Data(),fConversionCuts->GetEtaShift());
 221          fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 222       }
 223    }
 224
 225    return kTRUE;
 226 }
 227 //________________________________________________________________________
 228 void AliV0ReaderV1::UserExec(Option_t *option){
 229
 230   AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
 231   if(esdEvent) {
 232     if (!TGeoGlobalMagField::Instance()->GetField()) esdEvent->InitMagneticField();
 233   }
 234
 235
 236
 237    // Check if correctly initialized
 238    if(!fConversionGammas)Init();
 239
 240    // User Exec
 241    fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
 242
 243    // AOD Output
 244    FillAODOutput();
 245
 246 }
 247
 248 //________________________________________________________________________
 249 Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
 250 {
 251
 252    //Reset the TClonesArray
 253    fConversionGammas->Delete();
 254
 255    fInputEvent=inputEvent;
 256    fMCEvent=mcEvent;
 257
 258    if(!fInputEvent){
 259       AliError("No Input event");
 260       return kFALSE;
 261    }
 262    if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
 263
 264    // Count Primary Tracks Event
 265    CountTracks();
 266
 267    // Event Cuts
 268    if(!fConversionCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
 269
 270    // Set Magnetic Field
 271    AliKFParticle::SetField(fInputEvent->GetMagneticField());
 272
 273    if(fInputEvent->IsA()==AliESDEvent::Class()){
 274       ProcessESDV0s();
 275    }
 276    if(fInputEvent->IsA()==AliAODEvent::Class()){
 277       GetAODConversionGammas();
 278    }
 279
 280
 281
 282    return kTRUE;
 283 }
 284 ///________________________________________________________________________
 285 void AliV0ReaderV1::FillAODOutput()
 286 {
 287    // Fill AOD Output with reconstructed Photons
 288
 289    if(fInputEvent->IsA()==AliESDEvent::Class()){
 290       ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
 291       if(fCreateAOD) {
 292          AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 293          if (aodhandler && aodhandler->GetFillAOD()) {
 294             AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
 295             //PostData(0, fConversionGammas);
 296
 297          }
 298       }
 299    }
 300 }
 301
 302 ///________________________________________________________________________
 303 const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
 304
 305    // Get External Track Parameter with given charge
 306
 307    if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
 308
 309    // Check for sign flip
 310    if(fCurrentV0){
 311       if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
 312       if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
 313       if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
 314          tracklabel=fCurrentV0->GetPindex();
 315          return fCurrentV0->GetParamP();}
 316       if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
 317          tracklabel=fCurrentV0->GetNindex();
 318          return fCurrentV0->GetParamN();}
 319    }
 320    return 0x0;
 321 }
 322
 323 ///________________________________________________________________________
 324 Bool_t AliV0ReaderV1::ProcessESDV0s()
 325 {
 326    // Process ESD V0s for conversion photon reconstruction
 327    AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
 328
 329    AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
 330
 331    if(fESDEvent){
 332
 333       for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
 334          AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
 335          if(!fCurrentV0){
 336             printf("Requested V0 does not exist");
 337             continue;
 338          }
 339
 340          fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
 341
 342          if(fCurrentMotherKFCandidate){
 343             // Add Gamma to the TClonesArray
 344
 345             if(kUseAODConversionPhoton){
 346                new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
 347                AliAODConversionPhoton * currentConversionPhoton = (AliAODConversionPhoton*)(fConversionGammas->At(fConversionGammas->GetEntriesFast()-1));
 348                currentConversionPhoton->SetMass(fCurrentMotherKFCandidate->M());
 349                currentConversionPhoton->SetMassToZero();
 350
 351             }
 352             else{
 353                new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
 354             }
 355
 356             delete fCurrentMotherKFCandidate;
 357             fCurrentMotherKFCandidate=NULL;
 358          }
 359       }
 360    }
 361    return kTRUE;
 362 }
 363
 364 ///________________________________________________________________________
 365 AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
 366 {
 367         //   cout << currentV0Index << endl;
 368         // Reconstruct conversion photon from ESD v0
 369         fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonIn);
 370
 371         //checks if on the fly mode is set
 372         if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
 373                 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kOnFly);
 374                 return 0x0;
 375         }
 376
 377         // TrackLabels
 378         Int_t currentTrackLabels[2]={-1,-1};
 379
 380         // Get Daughter KF Particles
 381
 382         const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
 383         //    cout << fCurrentExternalTrackParamPositive << "\t" << currentTrackLabels[0] << endl;
 384         const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
 385         //    cout << fCurrentExternalTrackParamNegative << "\t" << currentTrackLabels[1] << endl;
 386         if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
 387
 388         // Apply some Cuts before Reconstruction
 389
 390         AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
 391         AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
 392         if(!negTrack || !posTrack) {
 393                 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kNoTracks);
 394                 return 0x0;
 395         }
 396         // Track Cuts
 397         if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
 398                 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kTrackCuts);
 399                 return 0x0;
 400         }
 401
 402         fConversionCuts->FillV0EtaBeforedEdxCuts(fCurrentV0->Eta());
 403         if (!fConversionCuts->dEdxCuts(posTrack)) {
 404                 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
 405                 return 0x0;
 406         }
 407         // PID Cuts
 408         if(!fConversionCuts->dEdxCuts(negTrack)) {
 409                 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
 410                 return 0x0;
 411         }
 412         fConversionCuts->FillV0EtaAfterdEdxCuts(fCurrentV0->Eta());
 413         // Reconstruct Photon
 414         AliKFConversionPhoton *fCurrentMotherKF=NULL;
 415         //    fUseConstructGamma = kFALSE;
 416         //    cout << "construct gamma " << endl;
 417         AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
 418         //    cout << fCurrentExternalTrackParamNegative << "\t" << endl;
 419         AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
 420         //    cout << fCurrentExternalTrackParamPositive << "\t"  << endl;
 421         //    cout << currentTrackLabels[0] << "\t" << currentTrackLabels[1] << endl;
 422         //    cout << "construct gamma " <<fUseConstructGamma << endl;
 423
 424         // Reconstruct Gamma
 425         if(fUseConstructGamma){
 426                 fCurrentMotherKF = new AliKFConversionPhoton();
 427                 fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
 428         }else{
 429                 fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
 430                 fCurrentMotherKF->SetMassConstraint(0,0.0001);
 431         }
 432
 433         // Set Track Labels
 434
 435         fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
 436
 437         // Set V0 index
 438
 439         fCurrentMotherKF->SetV0Index(currentV0Index);
 440
 441         //Set MC Label
 442         if(fMCEvent){
 443
 444                 AliStack *fMCStack= fMCEvent->Stack();
 445
 446                 Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
 447                 Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
 448
 449                 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
 450                 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
 451
 452                 if(fPositiveMCParticle&&fNegativeMCParticle){
 453                         fCurrentMotherKF->SetMCLabelPositive(labelp);
 454                         fCurrentMotherKF->SetMCLabelNegative(labeln);
 455                 }
 456         }
 457
 458         // Update Vertex (moved for same eta compared to old)
 459         //      cout << currentV0Index <<" \t before: \t" << fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz()  << endl;
 460         if(fUseImprovedVertex == kTRUE){
 461                 AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
 462                 //        cout << "Prim Vtx: " << primaryVertexImproved.GetX() << "\t" << primaryVertexImproved.GetY() << "\t" << primaryVertexImproved.GetZ() << endl;
 463                 primaryVertexImproved+=*fCurrentMotherKF;
 464                 fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
 465         }
 466         // SetPsiPair
 467
 468         Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
 469         fCurrentMotherKF->SetPsiPair(PsiPair);
 470
 471         // Recalculate ConversionPoint
 472         Double_t dca[2]={0,0};
 473         if(fUseOwnXYZCalculation){
 474                 Double_t convpos[3]={0,0,0};
 475                 if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
 476                         fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kConvPointFail);
 477                         delete fCurrentMotherKF;
 478                         fCurrentMotherKF=NULL;
 479                         return 0x0;
 480                 }
 481
 482                 fCurrentMotherKF->SetConversionPoint(convpos);
 483         }
 484
 485         if(fCurrentMotherKF->GetNDF() > 0.)
 486                 fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF());   //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
 487
 488
 489         // Set Dilepton Mass (moved down for same eta compared to old)
 490         fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
 491
 492         // Apply Photon Cuts
 493
 494         if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
 495                 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonCuts);
 496                 delete fCurrentMotherKF;
 497                 fCurrentMotherKF=NULL;
 498                 return 0x0;
 499         }
 500
 501         //    cout << currentV0Index <<" \t after: \t" <<fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz()  << endl;
 502
 503         fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonOut);
 504         return fCurrentMotherKF;
 505 }
 506
 507 ///________________________________________________________________________
 508 Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
 509    //
 510    // Angle between daughter momentum plane and plane
 511    //
 512
 513    AliExternalTrackParam nt(*negativeparam);
 514    AliExternalTrackParam pt(*positiveparam);
 515
 516    Float_t magField = fInputEvent->GetMagneticField();
 517
 518    Double_t xyz[3] = {0.,0.,0.};
 519    v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
 520
 521    // Double_t pPlus[3]  = {pt.Px(),pt.Py(),pt.Pz()};
 522    // Double_t pMinus[3] = {nt.Px(),nt.Py(),nt.Pz()};
 523
 524    // Double_t u[3] = {pPlus[0]+pMinus[0],pPlus[1]+pMinus[1],pPlus[2]+pMinus[2]};
 525    // Double_t normu = sqrt( (u[0]*u[0]) + (u[1]*u[1]) + (u[2]*u[2]) );
 526
 527    // u[0] = u[0] / normu;
 528    // u[1] = u[1] / normu;
 529    // u[2] = u[2] / normu;
 530
 531    // Double_t normpPlus = sqrt( (pPlus[0]*pPlus[0]) + (pPlus[1]*pPlus[1]) + (pPlus[2]*pPlus[2]) );
 532    // Double_t normpMinus = sqrt( (pMinus[0]*pMinus[0]) + (pMinus[1]*pMinus[1]) + (pMinus[2]*pMinus[2]) );
 533
 534    // pPlus[0] = pPlus[0] / normpPlus;
 535    // pPlus[1] = pPlus[1] / normpPlus;
 536    // pPlus[2] = pPlus[2] / normpPlus;
 537
 538    // pMinus[0] = pMinus[0] / normpMinus;
 539    // pMinus[1] = pMinus[1] / normpMinus;
 540    // pMinus[2] = pMinus[2] / normpMinus;
 541
 542    // Double_t v[3] = {0,0,0}; // pPlus X pMinus
 543    // v[0] = (pPlus[1]*pMinus[2]) - (pPlus[2]*pMinus[1]);
 544    // v[1] = (pPlus[2]*pMinus[0]) - (pPlus[0]*pMinus[2]);
 545    // v[2] = (pPlus[0]*pMinus[1]) - (pPlus[1]*pMinus[0]);
 546
 547    // Double_t w[3] = {0,0,0}; // u X v
 548    // w[0] = (u[1]*v[2]) - (u[2]*v[1]);
 549    // w[1] = (u[2]*v[0]) - (u[0]*v[2]);
 550    // w[2] = (u[0]*v[1]) - (u[1]*v[0]);
 551
 552    // Double_t z[3] = {0,0,1};
 553    // Double_t wc[3] = {0,0,0}; // u X v
 554    // wc[0] = (u[1]*z[2]) - (u[2]*z[1]);
 555    // wc[1] = (u[2]*z[0]) - (u[0]*z[2]);
 556    // wc[2] = (u[0]*z[1]) - (u[1]*z[0]);
 557
 558    // Double_t PhiV = TMath::ACos((w[0]*wc[0]) + (w[1]*wc[1]) + (w[2]*wc[2]));
 559    //return abs(PhiV);
 560
 561
 562    // TVector3 pPlus(pt.Px(),pt.Py(),pt.Pz());
 563    // TVector3 pMinus(nt.Px(),nt.Py(),nt.Pz());
 564
 565    // TVector3 u = pMinus + pPlus;
 566    // u = u*(1/u.Mag());
 567
 568    // TVector3 pHPlus = pPlus*(1/pPlus.Mag());
 569    // TVector3 pHMinus = pMinus*(1/pMinus.Mag());
 570
 571    // TVector3 v = pHPlus.Cross(pHMinus);
 572    // TVector3 w = u.Cross(v);
 573    // TVector3 z(0,0,1);
 574    // TVector3 wc = u.Cross(z);
 575
 576    // Double_t PhiV = w * wc;
 577
 578    Double_t mn[3] = {0,0,0};
 579    Double_t mp[3] = {0,0,0};
 580
 581    v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
 582    v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
 583
 584    Double_t deltat = 1.;
 585    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
 586    Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
 587
 588    Double_t momPosProp[3] = {0,0,0};
 589    Double_t momNegProp[3] = {0,0,0};
 590
 591    Double_t psiPair = 4.;
 592    if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
 593    if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
 594
 595    pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
 596    nt.GetPxPyPz(momNegProp);
 597
 598    Double_t pEle =
 599       TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
 600
 601    Double_t pPos =
 602       TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
 603
 604    Double_t scalarproduct =
 605       momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
 606
 607    Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
 608
 609 //    psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
 610    psiPair = TMath::ASin(deltat/chipair);
 611    return psiPair;
 612 }
 613
 614 ///________________________________________________________________________
 615 Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
 616
 617    // Get Center of the helix track parametrization
 618
 619    Int_t charge=track->Charge();
 620    Double_t b=fInputEvent->GetMagneticField();
 621
 622    Double_t     helix[6];
 623    track->GetHelixParameters(helix,b);
 624
 625    Double_t xpos =      helix[5];
 626    Double_t ypos =      helix[0];
 627    Double_t radius = TMath::Abs(1./helix[4]);
 628    Double_t phi = helix[2];
 629
 630    if(phi < 0){
 631       phi = phi + 2*TMath::Pi();
 632    }
 633
 634    phi -= TMath::Pi()/2.;
 635    Double_t xpoint =    radius * TMath::Cos(phi);
 636    Double_t ypoint =    radius * TMath::Sin(phi);
 637
 638    if(b<0){
 639       if(charge > 0){
 640          xpoint = - xpoint;
 641          ypoint = - ypoint;
 642       }
 643
 644       if(charge < 0){
 645          xpoint =       xpoint;
 646          ypoint =       ypoint;
 647       }
 648    }
 649    if(b>0){
 650       if(charge > 0){
 651          xpoint =       xpoint;
 652          ypoint =       ypoint;
 653       }
 654
 655       if(charge < 0){
 656          xpoint = - xpoint;
 657          ypoint = - ypoint;
 658       }
 659    }
 660    center[0] =  xpos + xpoint;
 661    center[1] =  ypos + ypoint;
 662
 663    return 1;
 664 }
 665 ///________________________________________________________________________
 666 Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
 667
 668    // Recalculate Conversion Point
 669
 670    if(!pparam||!nparam)return kFALSE;
 671
 672    Double_t helixcenterpos[2];
 673    GetHelixCenter(pparam,helixcenterpos);
 674
 675    Double_t helixcenterneg[2];
 676    GetHelixCenter(nparam,helixcenterneg);
 677
 678    Double_t helixpos[6];
 679    pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
 680    Double_t posradius = TMath::Abs(1./helixpos[4]);
 681
 682    Double_t helixneg[6];
 683    nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
 684    Double_t negradius = TMath::Abs(1./helixneg[4]);
 685
 686    // Calculate xy-position
 687
 688    Double_t xpos = helixcenterpos[0];
 689    Double_t ypos = helixcenterpos[1];
 690    Double_t xneg = helixcenterneg[0];
 691    Double_t yneg = helixcenterneg[1];
 692
 693    convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
 694    convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
 695
 696
 697    // Calculate Track XY vertex position
 698
 699    Double_t deltaXPos = convpos[0] -    xpos;
 700    Double_t deltaYPos = convpos[1] -    ypos;
 701
 702    Double_t deltaXNeg = convpos[0] -    xneg;
 703    Double_t deltaYNeg = convpos[1] -    yneg;
 704
 705    Double_t alphaPos =  TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
 706    Double_t alphaNeg =  TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
 707
 708    Double_t vertexXNeg =        xneg +  TMath::Abs(negradius)*TMath::Cos(alphaNeg);
 709    Double_t vertexYNeg =        yneg +  TMath::Abs(negradius)*TMath::Sin(alphaNeg);
 710
 711    Double_t vertexXPos =        xpos +  TMath::Abs(posradius)*TMath::Cos(alphaPos);
 712    Double_t vertexYPos =        ypos +  TMath::Abs(posradius)*TMath::Sin(alphaPos);
 713
 714    AliExternalTrackParam p(*pparam);
 715    AliExternalTrackParam n(*nparam);
 716
 717    TVector2 vertexPos(vertexXPos,vertexYPos);
 718    TVector2 vertexNeg(vertexXNeg,vertexYNeg);
 719
 720    // Convert to local coordinate system
 721    TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
 722    TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
 723
 724    // Propagate Track Params to Vertex
 725
 726    if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
 727    if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
 728
 729    // Check whether propagation was sucessful
 730
 731    if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
 732    if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
 733
 734    // Calculate z position
 735
 736    convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
 737
 738    // Calculate DCA
 739    TVector2 vdca=vertexPos-vertexNeg;
 740    dca[0]=vdca.Mod();
 741    dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
 742
 743    return kTRUE;
 744 }
 745 //________________________________________________________________________
 746 Bool_t AliV0ReaderV1::GetAODConversionGammas(){
 747
 748    // Get reconstructed conversion photons from satellite AOD file
 749
 750    AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
 751
 752    if(fAODEvent){
 753
 754       if(fConversionGammas == NULL){
 755          fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
 756       }
 757       fConversionGammas->Delete();//Reset the TClonesArray
 758
 759       //Get Gammas from satellite AOD gamma branch
 760
 761       AliAODConversionPhoton *gamma=0x0;
 762
 763       TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
 764
 765       if(!fInputGammas){
 766          FindDeltaAODBranchName();
 767          fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
 768       if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
 769       // Apply Selection Cuts to Gammas and create local working copy
 770       if(fInputGammas){
 771          for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
 772             gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
 773             if(gamma){
 774                if(fRelabelAODs)RelabelAODPhotonCandidates(gamma);
 775                if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
 776                   new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
 777             }
 778          }
 779       }
 780    }
 781
 782    if(fConversionGammas->GetEntries()){return kTRUE;}
 783
 784    return kFALSE;
 785 }
 786
 787 //________________________________________________________________________
 788 void AliV0ReaderV1::FindDeltaAODBranchName(){
 789
 790    // Find delta AOD branchname containing reconstructed photons
 791
 792    TList *list=fInputEvent->GetList();
 793    for(Int_t ii=0;ii<list->GetEntries();ii++){
 794       TString name((list->At(ii))->GetName());
 795       if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
 796          fDeltaAODBranchName=name;
 797          AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
 798          return;
 799       }
 800    }
 801 }
 802 //________________________________________________________________________
 803 void AliV0ReaderV1::RelabelAODPhotonCandidates(AliAODConversionPhoton *PhotonCandidate){
 804
 805    // Relabeling For AOD Event
 806    // ESDiD -> AODiD
 807    // MCLabel -> AODMCLabel
 808    Bool_t AODLabelPos = kFALSE;
 809    Bool_t AODLabelNeg = kFALSE;
 810
 811    for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
 812       AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
 813       if(!AODLabelPos){
 814          if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
 815             PhotonCandidate->SetMCLabelPositive(abs(tempDaughter->GetLabel()));
 816             PhotonCandidate->SetLabelPositive(i);
 817             AODLabelPos = kTRUE;
 818          }
 819       }
 820       if(!AODLabelNeg){
 821          if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
 822             PhotonCandidate->SetMCLabelNegative(abs(tempDaughter->GetLabel()));
 823             PhotonCandidate->SetLabelNegative(i);
 824             AODLabelNeg = kTRUE;
 825          }
 826       }
 827       if(AODLabelNeg && AODLabelPos){
 828          return;
 829       }
 830    }
 831    if(!AODLabelPos || !AODLabelNeg){
 832       AliError(Form("NO AOD Daughters Found Pos: %i %i Neg: %i %i",AODLabelPos,PhotonCandidate->GetTrackLabelPositive(),AODLabelNeg,PhotonCandidate->GetTrackLabelNegative()));
 833    }
 834
 835 }
 836
 837 //________________________________________________________________________
 838 void AliV0ReaderV1::CountTracks(){
 839
 840    // In Case of MC count only MB tracks
 841    // if(fMCEvent && fInputEvent->IsA()==AliESDEvent::Class()){
 842    //    fConversionCuts->GetNotRejectedParticles(1,NULL,fMCEvent);
 843    // }
 844    // if(fMCEvent && fInputEvent->IsA()==AliAODEvent::Class()){
 845    //    fConversionCuts->GetNotRejectedParticles(1,NULL,fInputEvent);
 846    // }
 847
 848    if(fInputEvent->IsA()==AliESDEvent::Class()){
 849       // Using standard function for setting Cuts
 850       Bool_t selectPrimaries=kTRUE;
 851       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 852       EsdTrackCuts->SetMaxDCAToVertexZ(2);
 853       EsdTrackCuts->SetEtaRange(-0.8, 0.8);
 854       EsdTrackCuts->SetPtRange(0.15);
 855       fNumberOfPrimaryTracks = 0;
 856       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 857          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 858          if(!curTrack) continue;
 859          if(!EsdTrackCuts->AcceptTrack(curTrack)) continue;
 860          //if(fMCEvent && !(fConversionCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()),fMCEvent->Stack(),fInputEvent))) continue;
 861          fNumberOfPrimaryTracks++;
 862       }
 863       delete EsdTrackCuts;
 864       EsdTrackCuts=0x0;
 865    }
 866    else if(fInputEvent->IsA()==AliAODEvent::Class()){
 867       fNumberOfPrimaryTracks = 0;
 868       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 869          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 870          if(curTrack->GetID()<0) continue; // Avoid double counting of tracks
 871          if(!curTrack->IsHybridGlobalConstrainedGlobal()) continue;
 872          if(abs(curTrack->Eta())>0.8) continue;
 873          if(curTrack->Pt()<0.15) continue;
 874          //if(fMCEvent && !(fConversionCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()),NULL,fInputEvent))) continue;
 875          //if(abs(curTrack->ZAtDCA())>2) continue; // Only Set For TPCOnly tracks
 876          fNumberOfPrimaryTracks++;
 877       }
 878    }
 879
 880    return;
 881 }
 882
 883 //________________________________________________________________________
 884 void AliV0ReaderV1::Terminate(Option_t *)
 885 {
 886
 887 }