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         fEventCuts(NULL),
  42         fConversionCuts(NULL),
  43         fTreeEvent(NULL),
  44         fTreeResolution(NULL),
  45         fPrimVtxZ(0.),
  46         fNContrVtx(0),
  47         fNESDtracksEta09(0),
  48         fNESDtracksEta0914(0),
  49         fNESDtracksEta14(0),
  50         fGammaRecCoords(5),
  51         fGammaMCCoords(5),
  52         fChi2ndf(0),
  53         fIsHeavyIon(0),
  54         fIsMC(kFALSE),
  55         fOutputList(NULL),
  56         fEventList(NULL),
  57         fResolutionList(NULL),
  58         fESDEvent(NULL),
  59         fMCEvent(NULL)
  60 {
  61
  62
  63 }
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskResolution::AliAnalysisTaskResolution(const char *name) : AliAnalysisTaskSE(name),
  67         fV0Reader(NULL),
  68         fConversionGammas(NULL),
  69         fEventCuts(NULL),
  70         fConversionCuts(NULL),
  71         fTreeEvent(NULL),
  72         fTreeResolution(NULL),
  73         fPrimVtxZ(0.),
  74         fNContrVtx(0),
  75         fNESDtracksEta09(0),
  76         fNESDtracksEta0914(0),
  77         fNESDtracksEta14(0),
  78         fGammaRecCoords(5),
  79         fGammaMCCoords(5),
  80         fChi2ndf(0),
  81         fIsHeavyIon(0),
  82         fIsMC(kFALSE),
  83         fOutputList(NULL),
  84         fEventList(NULL),
  85         fResolutionList(NULL),
  86         fESDEvent(NULL),
  87         fMCEvent(NULL)
  88 {
  89         // Default constructor
  90
  91         DefineInput(0, TChain::Class());
  92         DefineOutput(1, TList::Class());
  93 }
  94
  95 //________________________________________________________________________
  96 AliAnalysisTaskResolution::~AliAnalysisTaskResolution()
  97 {
  98         // default deconstructor
  99
 100 }
 101 //________________________________________________________________________
 102 void AliAnalysisTaskResolution::UserCreateOutputObjects()
 103 {
 104         // Create User Output Objects
 105
 106         if(fOutputList != NULL){
 107                 delete fOutputList;
 108                 fOutputList = NULL;
 109         }
 110         if(fOutputList == NULL){
 111                 fOutputList = new TList();
 112                 fOutputList->SetOwner(kTRUE);
 113         }
 114
 115         fEventList = new TList();
 116         fEventList->SetName("EventList");
 117         fEventList->SetOwner(kTRUE);
 118         fOutputList->Add(fEventList);
 119
 120         fTreeEvent = new TTree("Event","Event");
 121         fTreeEvent->Branch("primVtxZ",&fPrimVtxZ,"fPrimVtxZ/F");
 122         fTreeEvent->Branch("nContrVtx",&fNContrVtx,"fNContrVtx/I");
 123         fTreeEvent->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
 124         fTreeEvent->Branch("nGoodTracksEta14",&fNESDtracksEta14,"fNESDtracksEta14/I");
 125         fEventList->Add(fTreeEvent);
 126
 127         if (fIsMC){
 128                 fResolutionList = new TList();
 129                 fResolutionList->SetName("ResolutionList");
 130                 fResolutionList->SetOwner(kTRUE);
 131                 fOutputList->Add(fResolutionList);
 132
 133                 fTreeResolution = new TTree("Resolution","Resolution");
 134                 fTreeResolution->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
 135                 fTreeResolution->Branch("RecCoords",&fGammaRecCoords);
 136                 fTreeResolution->Branch("MCCoords",&fGammaMCCoords);
 137                 fTreeResolution->Branch("chi2ndf",&fChi2ndf,"fChi2ndf/F");
 138                 fResolutionList->Add(fTreeResolution);
 139         }
 140         PostData(1, fOutputList);
 141 }
 142
 143 //________________________________________________________________________
 144 void AliAnalysisTaskResolution::UserExec(Option_t *){
 145
 146         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 147
 148         Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
 149         if(eventQuality != 0){// Event Not Accepted
 150                 return;
 151         }
 152         fESDEvent = (AliESDEvent*) InputEvent();
 153         if (fESDEvent==NULL) return;
 154         if(MCEvent()){
 155                 fMCEvent = MCEvent();
 156         }
 157
 158         if(MCEvent()){
 159                 // Process MC Particle
 160                 if(fEventCuts->GetSignalRejection() != 0){
 161 //              if(fESDEvent->IsA()==AliESDEvent::Class()){
 162                         fEventCuts->GetNotRejectedParticles(    fEventCuts->GetSignalRejection(),
 163                                                                                                         fEventCuts->GetAcceptedHeader(),
 164                                                                                                         fMCEvent);
 165 //              }
 166 //              else if(fInputEvent->IsA()==AliAODEvent::Class()){
 167 //                      fEventCuts->GetNotRejectedParticles(    fEventCuts->GetSignalRejection(),
 168 //                                                                                                      fEventCuts->GetAcceptedHeader(),
 169 //                                                                                                      fInputEvent);
 170 //              }
 171                 }
 172         }
 173
 174
 175         if(fIsHeavyIon > 0 && !fEventCuts->IsCentralitySelected(fESDEvent)) return;
 176         fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
 177         fNESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
 178         fNESDtracksEta14 = fNESDtracksEta09 + fNESDtracksEta0914;
 179         if(fESDEvent){
 180                 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 181                         fNContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
 182                 } else {
 183                         fNContrVtx = 0;
 184                 }
 185 //              else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 186 //                      if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 187 //                              fNContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
 188 //                      } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 189 //                              fNContrVtx = 0;
 190 //                      }
 191 //              }
 192         }
 193         fPrimVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
 194
 195 //      if (fIsHeavyIon == 2){
 196 //              if (!(fNESDtracksEta09 > 20 && fNESDtracksEta09 < 80)) return;
 197 //      }
 198
 199
 200         if (fTreeEvent){
 201                 fTreeEvent->Fill();
 202         }
 203
 204         fConversionGammas=fV0Reader->GetReconstructedGammas();
 205         ProcessPhotons();
 206         PostData(1, fOutputList);
 207 }
 208
 209
 210 ///________________________________________________________________________
 211 void AliAnalysisTaskResolution::ProcessPhotons(){
 212
 213         // Fill Histograms for QA and MC
 214         for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
 215                 AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
 216                 if (gamma ==NULL) continue;
 217                 if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
 218                 fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
 219                 fGammaRecCoords(0) = gamma->GetPhotonPt();
 220                 fGammaRecCoords(1) = gamma->GetPhotonPhi();
 221                 fGammaRecCoords(2) = gamma->GetPhotonEta();
 222                 fGammaRecCoords(3) = gamma->GetConversionRadius();
 223                 fGammaRecCoords(4) = gamma->GetConversionZ();
 224                 fChi2ndf = gamma->GetChi2perNDF();
 225                 if(MCEvent()){
 226 //                      cout << "generating MC stack"<< endl;
 227                         AliStack *fMCStack = fMCEvent->Stack();
 228                         if (!fMCStack) continue;
 229                         TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
 230                         TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
 231 //                      cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
 232                         if(fMCStack && fEventCuts->GetSignalRejection() != 0){
 233                                 Int_t isPosFromMBHeader = fEventCuts->IsParticleFromBGEvent(gamma->GetMCLabelPositive(), fMCStack, fESDEvent);
 234                                 Int_t isNegFromMBHeader = fEventCuts->IsParticleFromBGEvent(gamma->GetMCLabelNegative(), fMCStack, fESDEvent);
 235                                 if( (isNegFromMBHeader < 1) || (isPosFromMBHeader < 1)) continue;
 236                         }
 237
 238                         if(posDaughter == NULL || negDaughter == NULL){
 239                                 continue;
 240                         } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
 241                                 continue;
 242                         } else {
 243 //                              cout << "same mother" << endl;
 244                                 Int_t pdgCodePos;
 245                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 246                                 Int_t pdgCodeNeg;
 247                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 248 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 249                                 Int_t pdgCode;
 250                                 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
 251 //                              cout << "PDG code: " << pdgCode << endl;
 252                                 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
 253                                         continue;
 254                                 else if ( !(pdgCodeNeg==pdgCodePos)){
 255                                         TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
 256                                         if(pdgCode == 111)
 257                                                 continue;
 258                                         else if (pdgCode == 221)
 259                                                 continue;
 260                                         else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
 261                                                 if(pdgCode == 22){
 262                                                         fGammaMCCoords(0) = truePhotonCanditate->Pt();
 263                                                         fGammaMCCoords(1) = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
 264                                                         fGammaMCCoords(2) = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
 265                                                         fGammaMCCoords(3) = gamma->GetNegativeMCDaughter(fMCStack)->R();
 266                                                         fGammaMCCoords(4) = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
 267
 268                                                         if (fTreeResolution){
 269                                                                 fTreeResolution->Fill();
 270                                                         }
 271                                                 }
 272                                         } else continue; //garbage
 273                                 } else continue; //garbage
 274                         }
 275                 }
 276         }
 277 }
 278
 279 //________________________________________________________________________
 280 Int_t AliAnalysisTaskResolution::CountTracks09(){
 281         Int_t fNumberOfESDTracks = 0;
 282         if(fInputEvent->IsA()==AliESDEvent::Class()){
 283         // Using standard function for setting Cuts
 284
 285                 AliStack *fMCStack = NULL;
 286                 if (MCEvent()){
 287                         fMCStack= fMCEvent->Stack();
 288                         if (!fMCStack) return 0;
 289                 }
 290
 291                 Bool_t selectPrimaries=kTRUE;
 292                 AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 293                 EsdTrackCuts->SetMaxDCAToVertexZ(2);
 294                 EsdTrackCuts->SetEtaRange(-0.9, 0.9);
 295                 EsdTrackCuts->SetPtRange(0.15);
 296
 297                 for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 298                         AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 299                         if(!curTrack) continue;
 300                         if(EsdTrackCuts->AcceptTrack(curTrack) ){
 301                                 if (fMCEvent){
 302                                         if(fMCStack && fEventCuts->GetSignalRejection() != 0){
 303                                                 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
 304                                                 if( (isFromMBHeader < 1) ) continue;
 305                                         }
 306                                 }
 307                                 fNumberOfESDTracks++;
 308                         }
 309                 }
 310                 delete EsdTrackCuts;
 311                 EsdTrackCuts=0x0;
 312         } else if(fInputEvent->IsA()==AliAODEvent::Class()){
 313                 for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 314                         AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 315                         if(!curTrack->IsPrimaryCandidate()) continue;
 316                         if(abs(curTrack->Eta())>0.9) continue;
 317                         if(curTrack->Pt()<0.15) continue;
 318                         if(abs(curTrack->ZAtDCA())>2) continue;
 319                         fNumberOfESDTracks++;
 320                 }
 321         }
 322
 323         return fNumberOfESDTracks;
 324 }
 325
 326 //________________________________________________________________________
 327 Int_t AliAnalysisTaskResolution::CountTracks0914(){
 328         Int_t fNumberOfESDTracks = 0;
 329         if(fInputEvent->IsA()==AliESDEvent::Class()){
 330                 // Using standard function for setting Cuts
 331
 332                 AliStack *fMCStack = NULL;
 333                 if (MCEvent()){
 334                         fMCStack= fMCEvent->Stack();
 335                         if (!fMCStack) return 0;
 336                 }
 337
 338                 AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 339                 EsdTrackCuts->SetMaxDCAToVertexZ(5);
 340                 EsdTrackCuts->SetEtaRange(0.9, 1.4);
 341                 EsdTrackCuts->SetPtRange(0.15);
 342
 343                 for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 344                         AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 345                         if(!curTrack) continue;
 346                         if(EsdTrackCuts->AcceptTrack(curTrack) ){
 347                                 if (fMCEvent){
 348                                         if(fMCStack && fEventCuts->GetSignalRejection() != 0){
 349                                                 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
 350                                                 if( (isFromMBHeader < 1) ) continue;
 351                                         }
 352                                 }
 353                                 fNumberOfESDTracks++;
 354                         }
 355                 }
 356                 EsdTrackCuts->SetEtaRange(-1.4, -0.9);
 357                 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 358                         AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 359                         if(!curTrack) continue;
 360                         if(EsdTrackCuts->AcceptTrack(curTrack) ){
 361                                 if (fMCEvent){
 362                                         if(fMCStack && fEventCuts->GetSignalRejection() != 0){
 363                                                 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
 364                                                 if( (isFromMBHeader < 1) ) continue;
 365                                         }
 366                                 }
 367                                 fNumberOfESDTracks++;
 368                         }
 369                 }
 370                 delete EsdTrackCuts;
 371                 EsdTrackCuts=0x0;
 372         } else if(fInputEvent->IsA()==AliAODEvent::Class()){
 373                 for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 374                         AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 375                         if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
 376                         if(curTrack->Pt()<0.15) continue;
 377                         if(abs(curTrack->ZAtDCA())>5) continue;
 378                         fNumberOfESDTracks++;
 379                 }
 380         }
 381
 382         return fNumberOfESDTracks;
 383 }
 384
 385 //________________________________________________________________________
 386 void AliAnalysisTaskResolution::Terminate(Option_t *)
 387 {
 388
 389 }