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