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    fV0Reader(NULL),
  40    fConversionGammas(NULL),
  41    fConversionCuts(NULL),
  42    fStreamMaterial(NULL),
  43    fStreamResolution(NULL),
  44    fIsHeavyIon(kFALSE),
  45    fOutputList(NULL),
  46    fESDEvent(NULL),
  47    fMCEvent(NULL)
  48 {
  49    // Default constructor
  50
  51    DefineInput(0, TChain::Class());
  52    DefineOutput(1, TList::Class());
  53 }
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
  57 {
  58    // default deconstructor
  59    if(fStreamMaterial){
  60       delete fStreamMaterial;
  61       fStreamMaterial = 0x0;
  62    }
  63    if(fStreamResolution){
  64       delete fStreamResolution;
  65       fStreamResolution = 0x0;
  66    }
  67 }
  68 //________________________________________________________________________
  69 void AliAnalysisTaskMaterial::UserCreateOutputObjects()
  70 {
  71    // Create User Output Objects
  72
  73    if(fOutputList != NULL){
  74       delete fOutputList;
  75       fOutputList = NULL;
  76    }
  77    if(fOutputList == NULL){
  78       fOutputList = new TList();
  79       fOutputList->SetOwner(kTRUE);
  80    }
  81
  82    // V0 Reader Cuts
  83    TString cutnumber = fConversionCuts->GetCutNumber();
  84
  85    fStreamMaterial = new TTreeSRedirector(Form("GammaConvV1_Material_%s.root",cutnumber.Data()),"recreate");
  86          fStreamResolution = new TTreeSRedirector(Form("GammaConvV1_Resolution_%s.root",cutnumber.Data()),"recreate");
  87    PostData(1, fOutputList);
  88 }
  89
  90 //________________________________________________________________________
  91 void AliAnalysisTaskMaterial::UserExec(Option_t *){
  92
  93    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
  94
  95    Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
  96         if(eventQuality != 0){// Event Not Accepted
  97       return;
  98    }
  99    fESDEvent = (AliESDEvent*) InputEvent();
 100    if (fESDEvent==NULL) return;
 101    if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
 102         Int_t nESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
 103         Int_t nESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
 104         Int_t nESDtracksEta14; // Estimate Event Multiplicity
 105   nESDtracksEta14= nESDtracksEta09 + nESDtracksEta0914;
 106         Int_t nContrVtx;
 107         if(fESDEvent){
 108                 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 109                         nContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
 110                 } else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 111                         if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 112                                 nContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
 113
 114                         } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 115                                 nContrVtx = 0;
 116                         }
 117                 }
 118    }
 119         Float_t primVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
 120
 121         if (fStreamMaterial){
 122                 (*fStreamMaterial)<<"Event"
 123                                                 << "primVtxZ=" << primVtxZ
 124                                                 << "nContrVtx=" << nContrVtx
 125                                                 << "nGoodTracksEta09=" << nESDtracksEta09
 126                                                 << "nGoodTracksEta0914=" << nESDtracksEta0914
 127                                                 << "nGoodTracksEta14=" << nESDtracksEta14
 128                                                 << "\n";
 129         }
 130
 131    fConversionGammas=fV0Reader->GetReconstructedGammas();
 132         if(MCEvent()){
 133                 fMCEvent = MCEvent();
 134         }
 135    ProcessPhotons();
 136         if(MCEvent()){
 137                 ProcessMCPhotons();
 138         }
 139    PostData(1, fOutputList);
 140 }
 141
 142 ///________________________________________________________________________
 143 void AliAnalysisTaskMaterial::FillMCTree(Int_t stackPos){
 144         AliStack *MCStack = fMCEvent->Stack();
 145         TParticle* candidate = (TParticle *)MCStack->Particle(stackPos);
 146         if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kFALSE)){
 147                 Float_t gammaPt = candidate->Pt();
 148                 Float_t gammaTheta = candidate->Theta();
 149                 if (fStreamMaterial){
 150                         (*fStreamMaterial)<<"AllGamma"
 151                                                         << "pt=" << gammaPt
 152 //                                                      << "p=" << gammaP
 153                                                         << "theta=" << gammaTheta
 154                                                         << "\n";
 155                 }
 156         }
 157         if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
 158                 Float_t gammaPt = candidate->Pt();
 159                 Float_t gammaTheta = candidate->Theta();
 160                 TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter());
 161                 TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter());
 162                 TVectorF coord(5);
 163                 coord(0) = daughter1->Vx();
 164                 coord(1) = daughter1->Vy();
 165                 coord(2) = daughter1->Vz();
 166                 coord(3) = daughter1->R();
 167                 coord(4) = candidate->Phi();
 168                 TVectorF daughterProp(4);
 169                 if (daughter1-> GetPdgCode() < 0){
 170                         daughterProp(0) = daughter2->Pt();
 171                         daughterProp(1) = daughter2->Theta();
 172                         daughterProp(2) = daughter1->Pt();
 173                         daughterProp(3) = daughter1->Theta();
 174                 } else {
 175                         daughterProp(0) = daughter1->Pt();
 176                         daughterProp(1) = daughter1->Theta();
 177                         daughterProp(2) = daughter2->Pt();
 178                         daughterProp(3) = daughter2->Theta();
 179                 }
 180                 if (fStreamMaterial){
 181                         (*fStreamMaterial)<<"ConvGammaMC"
 182                                                         << "pt=" << gammaPt
 183 //                                                      << "p=" << gammaP
 184                                                         << "theta=" << gammaTheta
 185                                                         << "coord.=" << &coord
 186                                                         << "daughterProp.=" << &daughterProp
 187                                                         << "\n";
 188                 }
 189         } // Converted MC Gamma
 190 }
 191
 192 ///________________________________________________________________________
 193 void AliAnalysisTaskMaterial::ProcessMCPhotons(){
 194         // Loop over all primary MC particle
 195         AliStack *ffMCStack = fMCEvent->Stack();
 196         for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
 197                 TParticle* particle = (TParticle *)ffMCStack->Particle(i);
 198                 if (!particle) continue;
 199                 if (particle->GetPdgCode() == 111 && particle->GetFirstDaughter() >= ffMCStack->GetNprimary()){
 200 //                      cout << "Undecayed pi0 found with mother: " << particle->GetMother(0) << endl;
 201                         for (Int_t j = 0; j < 2 ; j++){
 202                                 FillMCTree(particle->GetDaughter(j));
 203                         }
 204                 } else {
 205                         FillMCTree(i);
 206                 }
 207         }
 208 }
 209
 210 ///________________________________________________________________________
 211 void AliAnalysisTaskMaterial::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 //              cout << "i=  " <<firstGammaIndex << " of "<< fConversionGammas->GetEntriesFast() << endl;
 219       Float_t gammaPt = gamma->GetPhotonPt();
 220 //              Float_t gammaP = gamma->GetPhotonP();
 221                 Float_t gammaTheta = gamma->GetPhotonTheta();
 222                 Float_t gammaEta = gamma->GetPhotonEta();
 223                 Float_t gammaChi2NDF = gamma->GetChi2perNDF();
 224 //              Float_t gammaX = gamma->GetConversionX();
 225 //              Float_t gammaY = gamma->GetConversionY();
 226                 Float_t gammaZ = gamma->GetConversionZ();
 227       Float_t gammaR = gamma->GetConversionRadius();
 228                 Float_t gammaPhi = gamma->GetPhotonPhi();
 229
 230                 TVectorF coord(5);
 231                 coord(0) = gamma->GetConversionX();
 232                 coord(1) = gamma->GetConversionY();
 233                 coord(2) = gamma->GetConversionZ();
 234                 coord(3) = gamma->GetConversionRadius();
 235                 coord(4) = gamma->GetPhotonPhi();
 236                 TVectorF daughterProp(4);
 237                 AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
 238       AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
 239                 daughterProp(0) = posTrack->Pt();
 240       daughterProp(1) = posTrack->Theta();
 241       daughterProp(2) = negTrack->Pt();
 242       daughterProp(3) = negTrack->Theta();
 243
 244                 UInt_t kind = 9;
 245                 if(fMCEvent){
 246 //                      cout << "generating MC stack"<< endl;
 247                         AliStack *fMCStack = fMCEvent->Stack();
 248                         if (!fMCStack) continue;
 249                         TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
 250                         TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
 251 //                      cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
 252
 253                         if(posDaughter == NULL || negDaughter == NULL){
 254                                 kind = 9; // garbage
 255 //                              cout << "one of the daughters not available" << endl;
 256                         } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
 257                                 // Not Same Mother == Combinatorial Bck
 258                                 kind = 1;
 259 //                              cout << "not the same mother" << endl;
 260                                 Int_t pdgCodePos;
 261                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 262                                 Int_t pdgCodeNeg;
 263                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 264 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 265                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
 266                                         kind = 10; //Electron Combinatorial
 267                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
 268                                         kind = 15; //direct Electron Combinatorial
 269
 270                                 if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
 271                                         kind = 11; //Pion Combinatorial
 272                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) ||
 273                                         (TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
 274                                         kind = 12; //Pion, Proton Combinatorics
 275                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) ||
 276                                         (TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
 277                                         kind = 13; //Pion, Electron Combinatorics
 278                                 if (TMath::Abs(pdgCodePos)==321 || TMath::Abs(pdgCodeNeg)==321)
 279                                         kind = 14; //Kaon combinatorics
 280
 281                         } else {
 282 //                              cout << "same mother" << endl;
 283                                 Int_t pdgCodePos;
 284                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 285                                 Int_t pdgCodeNeg;
 286                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 287 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 288                                 Int_t pdgCode;
 289                                 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
 290 //                              cout << "PDG code: " << pdgCode << endl;
 291                                 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
 292                                         kind = 2; // combinatorics from hadronic decays
 293                                 else if ( !(pdgCodeNeg==pdgCodePos)){
 294                                         TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
 295                                         Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
 296                                         if(pdgCode == 111)
 297                                                 kind = 3; // pi0 Dalitz
 298                                         else if (pdgCode == 221)
 299                                                 kind = 4; // eta Dalitz
 300                                         else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
 301                                                 if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
 302                                                         kind = 0; // primary photons
 303                                                 } else if (pdgCode == 22){
 304                                                         kind = 5; //secondary photons
 305                                                 }
 306                                                 if(pdgCode == 22){
 307                                                         Float_t mcPt   = truePhotonCanditate->Pt();
 308                                                         Float_t mcR = gamma->GetNegativeMCDaughter(fMCStack)->R();
 309                                                         Float_t mcZ = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
 310                                                         Float_t mcPhi = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
 311                                                         Float_t mcEta = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
 312                                                         if (fStreamResolution){
 313                                                                 (*fStreamResolution)<<"Resolution"
 314                                                                 << "ESDpt=" << gammaPt
 315                                                                 << "ESDphi=" << gammaPhi
 316                                                                 << "ESDeta=" << gammaEta
 317                                                                 << "ESDr="<< gammaR
 318                                                                 << "ESDz="<< gammaZ
 319                                                                 << "MCpt=" << mcPt
 320                                                                 << "MCphi=" << mcPhi
 321                                                                 << "MCeta=" << mcEta
 322                                                                 << "MCr="<< mcR
 323                                                                 << "MCz="<< mcZ
 324                                                                 << "chi2ndf=" << gammaChi2NDF
 325                                                                 << "\n";
 326                                                         }
 327                                                 }
 328                                         } else  kind = 9; //garbage
 329                                 } else kind = 9; //garbage
 330                         }
 331 //                      cout << gammaPt << "\t" << gammaP<< "\t" << gammaEta<<  "\t" <<gammaChi2NDF << "\t" << gammaX<<  "\t" <<gammaY << "\t" << gammaZ<< "\t" << gammaR<<  "\t" <<  gammaPhi<< "\t" <<kind << endl;
 332
 333                         if (fStreamMaterial){
 334                                 (*fStreamMaterial)<<"ConvPointRec"
 335                                                                 << "pt=" << gammaPt
 336                                                                 << "theta=" << gammaTheta
 337                                                                 << "chi2ndf=" << gammaChi2NDF
 338                                                                 << "kind=" << kind
 339                                                                 << "coord.=" << &coord
 340                                                                 << "daughterProp.=" << &daughterProp
 341                                                                 << "\n";
 342                         }
 343                 } else {
 344                                 if (fStreamMaterial){
 345                                 (*fStreamMaterial)<<"ConvPointRec"
 346                                                                 << "pt=" << gammaPt
 347                                                                 << "theta=" << gammaTheta
 348                                                                 << "chi2ndf=" << gammaChi2NDF
 349                                                                 << "coord.=" << &coord
 350                                                                 << "daughterProp.=" << &daughterProp
 351                                                                 << "\n";
 352                         }
 353                 }
 354         }
 355 }
 356
 357 //________________________________________________________________________
 358 Int_t AliAnalysisTaskMaterial::CountTracks09(){
 359    Int_t fNumberOfESDTracks = 0;
 360    if(fInputEvent->IsA()==AliESDEvent::Class()){
 361    // Using standard function for setting Cuts
 362       Bool_t selectPrimaries=kTRUE;
 363       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 364       EsdTrackCuts->SetMaxDCAToVertexZ(2);
 365       EsdTrackCuts->SetEtaRange(-0.9, 0.9);
 366       EsdTrackCuts->SetPtRange(0.15);
 367
 368       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 369          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 370          if(!curTrack) continue;
 371          // if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
 372          //    if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
 373          // }
 374          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 375       }
 376       delete EsdTrackCuts;
 377       EsdTrackCuts=0x0;
 378    }
 379    else if(fInputEvent->IsA()==AliAODEvent::Class()){
 380       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 381          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 382          if(!curTrack->IsPrimaryCandidate()) continue;
 383          if(abs(curTrack->Eta())>0.9) continue;
 384          if(curTrack->Pt()<0.15) continue;
 385          if(abs(curTrack->ZAtDCA())>2) continue;
 386          fNumberOfESDTracks++;
 387       }
 388    }
 389
 390    return fNumberOfESDTracks;
 391 }
 392
 393 Int_t AliAnalysisTaskMaterial::CountTracks0914(){
 394
 395    // Using standard function for setting Cuts ; We use TPCOnlyTracks for outer eta region
 396    //Bool_t selectPrimaries=kTRUE;
 397          //   EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 398     Int_t fNumberOfESDTracks = 0;
 399     if(fInputEvent->IsA()==AliESDEvent::Class()){
 400       // Using standard function for setting Cuts
 401       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 402       EsdTrackCuts->SetMaxDCAToVertexZ(5);
 403       EsdTrackCuts->SetEtaRange(0.9, 1.4);
 404       EsdTrackCuts->SetPtRange(0.15);
 405
 406       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 407          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 408          if(!curTrack) continue;
 409          // if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
 410          //    if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
 411          // }
 412          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 413       }
 414       EsdTrackCuts->SetEtaRange(-1.4, -0.9);
 415       for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 416          AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 417          if(!curTrack) continue;
 418          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 419       }
 420       delete EsdTrackCuts;
 421       EsdTrackCuts=0x0;
 422    }
 423    else if(fInputEvent->IsA()==AliAODEvent::Class()){
 424       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 425          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 426 //          if(!curTrack->IsPrimaryCandidate()) continue;
 427          if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
 428          if(curTrack->Pt()<0.15) continue;
 429          if(abs(curTrack->ZAtDCA())>5) continue;
 430          fNumberOfESDTracks++;
 431       }
 432    }
 433
 434    return fNumberOfESDTracks;
 435 }
 436
 437 //________________________________________________________________________
 438 void AliAnalysisTaskMaterial::Terminate(Option_t *)
 439 {
 440 //    if (fStreamMaterial){
 441 //       fStreamMaterial->GetFile()->Write();
 442 //    }
 443 //    if (fStreamResolution){
 444 //       fStreamResolution->GetFile()->Write();
 445 //    }
 446 }