PWGGA/GammaConv/AliAnalysisTaskMaterial.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 "AliAnalysisTaskMaterial.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(AliAnalysisTaskMaterial)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskMaterial::AliAnalysisTaskMaterial(const char *name) : AliAnalysisTaskSE(name),
  39    fConversionGammas(NULL),
  40    fConversionCuts(NULL),
  41    fStreamMaterial(NULL),
  42    fStreamResolution(NULL),
  43    fIsHeavyIon(kFALSE),
  44    fOutputList(NULL),
  45    fESDEvent(NULL),
  46    fMCEvent(NULL)
  47 {
  48    // Default constructor
  49
  50    DefineInput(0, TChain::Class());
  51    DefineOutput(1, TList::Class());
  52 }
  53
  54 //________________________________________________________________________
  55 AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
  56 {
  57    // default deconstructor
  58    if(fStreamMaterial){
  59       delete fStreamMaterial;
  60       fStreamMaterial = 0x0;
  61    }
  62    if(fStreamResolution){
  63       delete fStreamResolution;
  64       fStreamResolution = 0x0;
  65    }
  66 }
  67 //________________________________________________________________________
  68 void AliAnalysisTaskMaterial::UserCreateOutputObjects()
  69 {
  70    // Create User Output Objects
  71
  72    if(fOutputList != NULL){
  73       delete fOutputList;
  74       fOutputList = NULL;
  75    }
  76    if(fOutputList == NULL){
  77       fOutputList = new TList();
  78       fOutputList->SetOwner(kTRUE);
  79    }
  80
  81    // V0 Reader Cuts
  82    TString cutnumber = fConversionCuts->GetCutNumber();
  83
  84    fStreamMaterial = new TTreeSRedirector(Form("GammaConvV1_Material_%s.root",cutnumber.Data()));
  85         fStreamResolution = new TTreeSRedirector(Form("GammaConvV1_Resolution_%s.root",cutnumber.Data()));
  86    PostData(1, fOutputList);
  87 }
  88
  89 //________________________________________________________________________
  90 void AliAnalysisTaskMaterial::UserExec(Option_t *){
  91
  92    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
  93
  94    Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
  95    if(eventQuality != 0){// Event Not Accepted
  96       return;
  97    }
  98    fESDEvent = (AliESDEvent*) InputEvent();
  99    if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
 100         Int_t nESDtracks = CountESDTracks(); // Estimate Event Multiplicity
 101         Int_t nContrVtx;
 102         if(fESDEvent){
 103                 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 104                         nContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
 105                 } else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 106                         if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 107                                 nContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
 108
 109                         } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 110                                 nContrVtx = 0;
 111                         }
 112                 }
 113    }
 114         Int_t primVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
 115
 116         if (fStreamMaterial){
 117                 (*fStreamMaterial)<<"Event"
 118                                                 << "primVtxZ=" << primVtxZ
 119                                                 << "nContrVtx=" << nContrVtx
 120                                                 << "nGoodTracks=" << nESDtracks
 121                                                 << "\n";
 122         }
 123
 124    fConversionGammas=fV0Reader->GetReconstructedGammas();
 125         if(MCEvent()){
 126                 fMCEvent = MCEvent();
 127         }
 128    ProcessPhotons();
 129         if(MCEvent()){
 130                 ProcessMCPhotons();
 131         }
 132    PostData(1, fOutputList);
 133 }
 134
 135
 136 ///________________________________________________________________________
 137 void AliAnalysisTaskMaterial::ProcessMCPhotons(){
 138         // Loop over all primary MC particle
 139         AliStack *ffMCStack = fMCEvent->Stack();
 140         for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
 141                 TParticle* particle = (TParticle *)ffMCStack->Particle(i);
 142                 if (!particle) continue;
 143                 if(fConversionCuts->PhotonIsSelectedMC(particle,ffMCStack,kFALSE)){
 144                         Float_t gammaPt = particle->Pt();
 145                         Float_t gammaP = particle->P();
 146                         Float_t gammaEta = particle->Eta();
 147                         if (fStreamMaterial){
 148                                 (*fStreamMaterial)<<"AllGamma"
 149                                                                 << "pt=" << gammaPt
 150                                                                 << "p=" << gammaP
 151                                                                 << "eta=" << gammaEta
 152                                                                 << "\n";
 153                         }
 154                 }
 155                 if(fConversionCuts->PhotonIsSelectedMC(particle,ffMCStack,kTRUE)){
 156                         Float_t gammaPt = particle->Pt();
 157                         Float_t gammaP = particle->P();
 158                         Float_t gammaEta = particle->Eta();
 159                         TParticle* daughter = (TParticle *)ffMCStack->Particle(particle->GetFirstDaughter());
 160                         Float_t gammaX = daughter->Vx();
 161                         Float_t gammaY =  daughter->Vy();
 162                         Float_t gammaZ = daughter->Vz();
 163                         Float_t gammaR = daughter->R();
 164                         Float_t gammaPhi = particle->Phi();
 165
 166                         if (fStreamMaterial){
 167                                 (*fStreamMaterial)<<"ConvGammaMC"
 168                                                                 << "pt=" << gammaPt
 169                                                                 << "p=" << gammaP
 170                                                                 << "eta=" << gammaEta
 171                                                                 << "X=" << gammaX
 172                                                                 << "Y=" << gammaY
 173                                                                 << "Z=" << gammaZ
 174                                                                 << "R=" << gammaR
 175                                                                 << "Phi=" << gammaPhi
 176                                                                 << "\n";
 177                         }
 178                 } // Converted MC Gamma
 179         }
 180 }
 181
 182 ///________________________________________________________________________
 183 void AliAnalysisTaskMaterial::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(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
 189 //              cout << "i=  " <<firstGammaIndex << " of "<< fConversionGammas->GetEntriesFast() << endl;
 190       Float_t gammaPt = gamma->GetPhotonPt();
 191                 Float_t gammaP = gamma->GetPhotonP();
 192                 Float_t gammaEta = gamma->GetPhotonEta();
 193                 Float_t gammaChi2NDF = gamma->GetChi2perNDF();
 194                 Float_t gammaX = gamma->GetConversionX();
 195                 Float_t gammaY = gamma->GetConversionY();
 196                 Float_t gammaZ = gamma->GetConversionZ();
 197       Float_t gammaR = gamma->GetConversionRadius();
 198                 Float_t gammaPhi = gamma->GetPhotonPhi();
 199                 UInt_t kind = 9;
 200                 if(fMCEvent){
 201 //                      cout << "generating MC stack"<< endl;
 202                         AliStack *fMCStack = fMCEvent->Stack();
 203                         if (!fMCStack) continue;
 204                         TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
 205                         TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
 206 //                      cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
 207
 208                         if(posDaughter == NULL || negDaughter == NULL){
 209                                 kind = 9; // garbage
 210 //                              cout << "one of the daughters not available" << endl;
 211                         } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
 212                                 // Not Same Mother == Combinatorial Bck
 213                                 kind = 1;
 214 //                              cout << "not the same mother" << endl;
 215                                 Int_t pdgCodePos;
 216                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 217                                 Int_t pdgCodeNeg;
 218                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 219 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 220                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
 221                                         kind = 10; //Electron Combinatorial
 222                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
 223                                         kind = 15; //direct Electron Combinatorial
 224
 225                                 if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
 226                                         kind = 11; //Pion Combinatorial
 227                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) ||
 228                                         (TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
 229                                         kind = 12; //Pion, Proton Combinatorics
 230                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) ||
 231                                         (TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
 232                                         kind = 13; //Pion, Electron Combinatorics
 233                                 if (TMath::Abs(pdgCodePos)==321 || TMath::Abs(pdgCodeNeg)==321)
 234                                         kind = 14; //Kaon combinatorics
 235
 236                         } else {
 237 //                              cout << "same mother" << endl;
 238                                 Int_t pdgCodePos;
 239                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 240                                 Int_t pdgCodeNeg;
 241                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 242 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 243                                 Int_t pdgCode;
 244                                 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
 245 //                              cout << "PDG code: " << pdgCode << endl;
 246                                 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
 247                                         kind = 2; // combinatorics from hadronic decays
 248                                 else if ( !(pdgCodeNeg==pdgCodePos)){
 249                                         TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
 250                                         if(pdgCode == 111)
 251                                                 kind = 3; // pi0 Dalitz
 252                                         else if (pdgCode == 221)
 253                                                 kind = 4; // eta Dalitz
 254                                         else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
 255                                                 if(pdgCode == 22 && negDaughter->GetMother(0) <= fMCStack->GetNprimary()){
 256                                                         kind = 0; // primary photons
 257                                                 } else if (pdgCode == 22){
 258                                                         kind = 5; //secondary photons
 259                                                 }
 260                                                 if(pdgCode == 22){
 261                                                         Float_t mcPt   = truePhotonCanditate->Pt();
 262                                                         Float_t mcR = gamma->GetNegativeMCDaughter(fMCStack)->R();
 263                                                         Float_t mcZ = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
 264                                                         Float_t mcPhi = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
 265                                                         Float_t mcEta = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
 266                                                         if (fStreamResolution){
 267                                                                 (*fStreamResolution)<<"Resolution"
 268                                                                 << "ESDpt=" << gammaPt
 269                                                                 << "ESDphi=" << gammaPhi
 270                                                                 << "ESDeta=" << gammaEta
 271                                                                 << "ESDr="<< gammaR
 272                                                                 << "ESDz="<< gammaZ
 273                                                                 << "MCpt=" << mcPt
 274                                                                 << "MCphi=" << mcPhi
 275                                                                 << "MCeta=" << mcEta
 276                                                                 << "MCr="<< mcR
 277                                                                 << "MCz="<< mcZ
 278                                                                 << "chi2ndf=" << gammaChi2NDF
 279                                                                 << "\n";
 280                                                         }
 281                                                 }
 282                                         } else  kind = 9; //garbage
 283                                 } else kind = 9; //garbage
 284                         }
 285 //                      cout << gammaPt << "\t" << gammaP<< "\t" << gammaEta<<  "\t" <<gammaChi2NDF << "\t" << gammaX<<  "\t" <<gammaY << "\t" << gammaZ<< "\t" << gammaR<<  "\t" <<  gammaPhi<< "\t" <<kind << endl;
 286
 287                         if (fStreamMaterial){
 288                                 (*fStreamMaterial)<<"ConvPointRec"
 289                                                                 << "pt=" << gammaPt
 290                                                                 << "p=" << gammaP
 291                                                                 << "eta=" << gammaEta
 292                                                                 << "chi2ndf=" << gammaChi2NDF
 293                                                                 << "kind=" << kind
 294                                                                 << "X=" << gammaX
 295                                                                 << "Y=" << gammaY
 296                                                                 << "Z=" << gammaZ
 297                                                                 << "R=" << gammaR
 298                                                                 << "Phi=" << gammaPhi
 299                                                                 << "\n";
 300                         }
 301                 } else {
 302                                 if (fStreamMaterial){
 303                                 (*fStreamMaterial)<<"ConvPointRec"
 304                                                                 << "pt=" << gammaPt
 305                                                                 << "p=" << gammaP
 306                                                                 << "eta=" << gammaEta
 307                                                                 << "chi2ndf=" << gammaChi2NDF
 308                                                                 << "X=" << gammaX
 309                                                                 << "Y=" << gammaY
 310                                                                 << "Z=" << gammaZ
 311                                                                 << "R=" << gammaR
 312                                                                 << "Phi=" << gammaPhi
 313                                                                 << "\n";
 314                         }
 315                 }
 316         }
 317 }
 318
 319 //________________________________________________________________________
 320 Int_t AliAnalysisTaskMaterial::CountESDTracks(){
 321
 322    AliESDtrackCuts *EsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
 323    // Using standard function for setting Cuts
 324    Bool_t selectPrimaries=kTRUE;
 325    EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 326    EsdTrackCuts->SetMaxDCAToVertexZ(2);
 327    EsdTrackCuts->SetEtaRange(-0.8, 0.8);
 328    EsdTrackCuts->SetPtRange(0.15);
 329
 330    Int_t fNumberOfESDTracks = 0;
 331    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 332       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 333       if(!curTrack) continue;
 334       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 335    }
 336    delete EsdTrackCuts;
 337    EsdTrackCuts=0x0;
 338
 339    return fNumberOfESDTracks;
 340 }
 341
 342
 343 //________________________________________________________________________
 344 void AliAnalysisTaskMaterial::Terminate(Option_t *)
 345 {
 346    if (fStreamMaterial){
 347       fStreamMaterial->GetFile()->Write();
 348    }
 349    if (fStreamResolution){
 350       fStreamResolution->GetFile()->Write();
 351    }
 352 }