PWGGA/GammaConv/AliAnalysisTaskResolution.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                                                                                                                                              *
   4  * Authors: Friederike Bock                                                                                                                     *
   5  * Version 1.0                                                                                                                                                          *
   6  *                                                                                                                                                                                              *
   7  * Permission to use, copy, modify and distribute this software and its  *
   8  * documentation strictly for non-commercial purposes is hereby granted  *
   9  * without fee, provided that the above copyright notice appears in all  *
  10  * copies and that both the copyright notice and this permission notice  *
  11  * appear in the supporting documentation. The authors make no claims    *
  12  * about the suitability of this software for any purpose. It is                        *
  13  * provided "as is" without express or implied warranty.                                                *
  14  **************************************************************************/
  15
  16 ////////////////////////////////////////////////
  17 //---------------------------------------------
  18 // QA Task for V0 Reader V1
  19 //---------------------------------------------
  20 ////////////////////////////////////////////////
  21
  22 #include "AliAnalysisTaskResolution.h"
  23 #include "TChain.h"
  24 #include "AliAnalysisManager.h"
  25 #include "TParticle.h"
  26 #include "TVectorF.h"
  27 #include "AliPIDResponse.h"
  28 #include "AliESDtrackCuts.h"
  29 #include "TFile.h"
  30
  31 class iostream;
  32
  33 using namespace std;
  34
  35 ClassImp(AliAnalysisTaskResolution)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskResolution::AliAnalysisTaskResolution() : AliAnalysisTaskSE(),
  39    fV0Reader(NULL),
  40    fConversionGammas(NULL),
  41    fConversionCuts(NULL),
  42    fTreeEvent(NULL),
  43    fTreeResolution(NULL),
  44    fPrimVtxZ(0.),
  45    fNContrVtx(0),
  46    fNESDtracksEta09(0),
  47    fNESDtracksEta0914(0),
  48    fNESDtracksEta14(0),
  49    fGammaRecCoords(5),
  50    fGammaMCCoords(5),
  51    fChi2ndf(0),
  52    fIsHeavyIon(kFALSE),
  53    fOutputList(NULL),
  54    fEventList(NULL),
  55    fResolutionList(NULL),
  56    fESDEvent(NULL),
  57    fMCEvent(NULL)
  58 {
  59
  60
  61 }
  62
  63
  64 //________________________________________________________________________
  65 AliAnalysisTaskResolution::AliAnalysisTaskResolution(const char *name) : AliAnalysisTaskSE(name),
  66    fV0Reader(NULL),
  67    fConversionGammas(NULL),
  68    fConversionCuts(NULL),
  69    fTreeEvent(NULL),
  70    fTreeResolution(NULL),
  71    fPrimVtxZ(0.),
  72    fNContrVtx(0),
  73    fNESDtracksEta09(0),
  74    fNESDtracksEta0914(0),
  75    fNESDtracksEta14(0),
  76    fGammaRecCoords(5),
  77    fGammaMCCoords(5),
  78    fChi2ndf(0),
  79    fIsHeavyIon(kFALSE),
  80    fOutputList(NULL),
  81    fEventList(NULL),
  82    fResolutionList(NULL),
  83    fESDEvent(NULL),
  84    fMCEvent(NULL)
  85 {
  86    // Default constructor
  87
  88    DefineInput(0, TChain::Class());
  89    DefineOutput(1, TList::Class());
  90 }
  91
  92 //________________________________________________________________________
  93 AliAnalysisTaskResolution::~AliAnalysisTaskResolution()
  94 {
  95    // default deconstructor
  96
  97 }
  98 //________________________________________________________________________
  99 void AliAnalysisTaskResolution::UserCreateOutputObjects()
 100 {
 101    // Create User Output Objects
 102
 103    if(fOutputList != NULL){
 104       delete fOutputList;
 105       fOutputList = NULL;
 106    }
 107    if(fOutputList == NULL){
 108       fOutputList = new TList();
 109       fOutputList->SetOwner(kTRUE);
 110    }
 111
 112    fEventList = new TList();
 113    fEventList->SetName("EventList");
 114    fEventList->SetOwner(kTRUE);
 115    fOutputList->Add(fEventList);
 116
 117    fTreeEvent = new TTree("Event","Event");
 118    fTreeEvent->Branch("primVtxZ",&fPrimVtxZ,"fPrimVtxZ/F");
 119    fTreeEvent->Branch("nContrVtx",&fNContrVtx,"fNContrVtx/I");
 120    fTreeEvent->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
 121    fTreeEvent->Branch("nGoodTracksEta0914",&fNESDtracksEta0914,"fNESDtracksEta0914/I");
 122    fTreeEvent->Branch("nGoodTracksEta14",&fNESDtracksEta14,"fNESDtracksEta14/I");
 123    fEventList->Add(fTreeEvent);
 124
 125    fResolutionList = new TList();
 126    fResolutionList->SetName("ResolutionList");
 127    fResolutionList->SetOwner(kTRUE);
 128    fOutputList->Add(fResolutionList);
 129
 130    fTreeResolution = new TTree("Resolution","Resolution");
 131    fTreeResolution->Branch("RecCoords",&fGammaRecCoords);
 132    fTreeResolution->Branch("MCCoords",&fGammaMCCoords);
 133    fTreeResolution->Branch("chi2ndf",&fChi2ndf,"fChi2ndf/F");
 134    fResolutionList->Add(fTreeResolution);
 135    // V0 Reader Cuts
 136    TString cutnumber = fConversionCuts->GetCutNumber();
 137
 138    PostData(1, fOutputList);
 139 }
 140
 141 //________________________________________________________________________
 142 void AliAnalysisTaskResolution::UserExec(Option_t *){
 143
 144    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 145
 146       Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
 147    if(eventQuality != 0){// Event Not Accepted
 148       return;
 149    }
 150    fESDEvent = (AliESDEvent*) InputEvent();
 151    if (fESDEvent==NULL) return;
 152    if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
 153    fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
 154    fNESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
 155    fNESDtracksEta14 = fNESDtracksEta09 + fNESDtracksEta0914;
 156    if(fESDEvent){
 157       if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 158          fNContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
 159       } else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 160          if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 161             fNContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
 162          } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 163             fNContrVtx = 0;
 164          }
 165       }
 166    }
 167    fPrimVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
 168
 169    if (fTreeEvent){
 170       fTreeEvent->Fill();
 171    }
 172
 173    fConversionGammas=fV0Reader->GetReconstructedGammas();
 174         if(MCEvent()){
 175                 fMCEvent = MCEvent();
 176         }
 177    ProcessPhotons();
 178    PostData(1, fOutputList);
 179 }
 180
 181
 182 ///________________________________________________________________________
 183 void AliAnalysisTaskResolution::ProcessPhotons(){
 184
 185         // Fill Histograms for QA and MC
 186    for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
 187       AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
 188       if (gamma ==NULL) continue;
 189       if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
 190       fGammaRecCoords(0) = gamma->GetPhotonPt();
 191       fGammaRecCoords(1) = gamma->GetPhotonPhi();
 192       fGammaRecCoords(2) = gamma->GetPhotonEta();
 193       fGammaRecCoords(3) = gamma->GetConversionRadius();
 194       fGammaRecCoords(4) = gamma->GetConversionZ();
 195       fChi2ndf = gamma->GetChi2perNDF();
 196                 if(fMCEvent){
 197 //                      cout << "generating MC stack"<< endl;
 198                         AliStack *fMCStack = fMCEvent->Stack();
 199                         if (!fMCStack) continue;
 200                         TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
 201                         TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
 202 //                      cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
 203
 204                         if(posDaughter == NULL || negDaughter == NULL){
 205                                 continue;
 206                         } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
 207                                 continue;
 208                         } else {
 209 //                              cout << "same mother" << endl;
 210                                 Int_t pdgCodePos;
 211                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 212                                 Int_t pdgCodeNeg;
 213                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 214 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 215                                 Int_t pdgCode;
 216                                 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
 217 //                              cout << "PDG code: " << pdgCode << endl;
 218                                 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
 219                                         continue;
 220                                 else if ( !(pdgCodeNeg==pdgCodePos)){
 221                                         TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
 222                                         if(pdgCode == 111)
 223                                                 continue;
 224                                         else if (pdgCode == 221)
 225                                                 continue;
 226                                         else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
 227                                                 if(pdgCode == 22){
 228                      fGammaMCCoords(0) = truePhotonCanditate->Pt();
 229                      fGammaMCCoords(1) = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
 230                      fGammaMCCoords(2) = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
 231                      fGammaMCCoords(3) = gamma->GetNegativeMCDaughter(fMCStack)->R();
 232                      fGammaMCCoords(4) = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
 233
 234                      if (fTreeResolution){
 235                         fTreeResolution->Fill();
 236                                                         }
 237                                                 }
 238                                         } else continue; //garbage
 239                                 } else continue; //garbage
 240                         }
 241       }
 242         }
 243 }
 244
 245 //________________________________________________________________________
 246 Int_t AliAnalysisTaskResolution::CountTracks09(){
 247    Int_t fNumberOfESDTracks = 0;
 248    if(fInputEvent->IsA()==AliESDEvent::Class()){
 249    // Using standard function for setting Cuts
 250       Bool_t selectPrimaries=kTRUE;
 251       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 252       EsdTrackCuts->SetMaxDCAToVertexZ(2);
 253       EsdTrackCuts->SetEtaRange(-0.9, 0.9);
 254       EsdTrackCuts->SetPtRange(0.15);
 255
 256       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 257          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 258          if(!curTrack) continue;
 259          // if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
 260          //    if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
 261          // }
 262          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 263       }
 264       delete EsdTrackCuts;
 265       EsdTrackCuts=0x0;
 266    }
 267    else if(fInputEvent->IsA()==AliAODEvent::Class()){
 268       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 269          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 270          if(!curTrack->IsPrimaryCandidate()) continue;
 271          if(abs(curTrack->Eta())>0.9) continue;
 272          if(curTrack->Pt()<0.15) continue;
 273          if(abs(curTrack->ZAtDCA())>2) continue;
 274          fNumberOfESDTracks++;
 275       }
 276    }
 277
 278    return fNumberOfESDTracks;
 279 }
 280
 281 Int_t AliAnalysisTaskResolution::CountTracks0914(){
 282
 283    // Using standard function for setting Cuts ; We use TPCOnlyTracks for outer eta region
 284    //Bool_t selectPrimaries=kTRUE;
 285          //   EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 286     Int_t fNumberOfESDTracks = 0;
 287     if(fInputEvent->IsA()==AliESDEvent::Class()){
 288       // Using standard function for setting Cuts
 289       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 290       EsdTrackCuts->SetMaxDCAToVertexZ(5);
 291       EsdTrackCuts->SetEtaRange(0.9, 1.4);
 292       EsdTrackCuts->SetPtRange(0.15);
 293
 294       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 295          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 296          if(!curTrack) continue;
 297          // if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
 298          //    if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
 299          // }
 300          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 301       }
 302       EsdTrackCuts->SetEtaRange(-1.4, -0.9);
 303       for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 304          AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 305          if(!curTrack) continue;
 306          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 307       }
 308       delete EsdTrackCuts;
 309       EsdTrackCuts=0x0;
 310    }
 311    else if(fInputEvent->IsA()==AliAODEvent::Class()){
 312       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 313          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 314 //          if(!curTrack->IsPrimaryCandidate()) continue;
 315          if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
 316          if(curTrack->Pt()<0.15) continue;
 317          if(abs(curTrack->ZAtDCA())>5) continue;
 318          fNumberOfESDTracks++;
 319       }
 320    }
 321
 322    return fNumberOfESDTracks;
 323 }
 324
 325 //________________________________________________________________________
 326 void AliAnalysisTaskResolution::Terminate(Option_t *)
 327 {
 328 //    if (fStreamMaterial){
 329 //       fStreamMaterial->GetFile()->Write();
 330 //    }
 331 //    if (fStreamResolution){
 332 //       fStreamResolution->GetFile()->Write();
 333 //    }
 334 }