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 = CountESDTracks09(); // Estimate Event Multiplicity
 103         Int_t nESDtracksEta0914 = CountESDTracks0914(); // Estimate Event Multiplicity
 104         //      Int_t nESDtracksEta14 = CountESDTracks14(); // Estimate Event Multiplicity
 105         Int_t nESDtracksEta14; // Estimate Event Multiplicity
 106   nESDtracksEta14= nESDtracksEta09 + nESDtracksEta0914;
 107         Int_t nContrVtx;
 108         if(fESDEvent){
 109                 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 110                         nContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
 111                 } else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 112                         if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 113                                 nContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
 114
 115                         } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 116                                 nContrVtx = 0;
 117                         }
 118                 }
 119    }
 120         Float_t primVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
 121
 122         if (fStreamMaterial){
 123                 (*fStreamMaterial)<<"Event"
 124                                                 << "primVtxZ=" << primVtxZ
 125                                                 << "nContrVtx=" << nContrVtx
 126                                                 << "nGoodTracksEta09=" << nESDtracksEta09
 127                                                 << "nGoodTracksEta0914=" << nESDtracksEta0914
 128                                                 << "nGoodTracksEta14=" << nESDtracksEta14
 129                                                 << "\n";
 130         }
 131
 132    fConversionGammas=fV0Reader->GetReconstructedGammas();
 133         if(MCEvent()){
 134                 fMCEvent = MCEvent();
 135         }
 136    ProcessPhotons();
 137         if(MCEvent()){
 138                 ProcessMCPhotons();
 139         }
 140    PostData(1, fOutputList);
 141 }
 142
 143 ///________________________________________________________________________
 144 void AliAnalysisTaskMaterial::FillMCTree(Int_t stackPos){
 145         AliStack *MCStack = fMCEvent->Stack();
 146         TParticle* candidate = (TParticle *)MCStack->Particle(stackPos);
 147         if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kFALSE)){
 148                 Float_t gammaPt = candidate->Pt();
 149                 Float_t gammaTheta = candidate->Theta();
 150                 if (fStreamMaterial){
 151                         (*fStreamMaterial)<<"AllGamma"
 152                                                         << "pt=" << gammaPt
 153 //                                                      << "p=" << gammaP
 154                                                         << "theta=" << gammaTheta
 155                                                         << "\n";
 156                 }
 157         }
 158         if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
 159                 Float_t gammaPt = candidate->Pt();
 160                 Float_t gammaTheta = candidate->Theta();
 161                 TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter());
 162                 TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter());
 163                 TVectorF coord(5);
 164                 coord(0) = daughter1->Vx();
 165                 coord(1) = daughter1->Vy();
 166                 coord(2) = daughter1->Vz();
 167                 coord(3) = daughter1->R();
 168                 coord(4) = candidate->Phi();
 169                 TVectorF daughterProp(4);
 170                 if (daughter1-> GetPdgCode() < 0){
 171                         daughterProp(0) = daughter2->Pt();
 172                         daughterProp(1) = daughter2->Theta();
 173                         daughterProp(2) = daughter1->Pt();
 174                         daughterProp(3) = daughter1->Theta();
 175                 } else {
 176                         daughterProp(0) = daughter1->Pt();
 177                         daughterProp(1) = daughter1->Theta();
 178                         daughterProp(2) = daughter2->Pt();
 179                         daughterProp(3) = daughter2->Theta();
 180                 }
 181                 if (fStreamMaterial){
 182                         (*fStreamMaterial)<<"ConvGammaMC"
 183                                                         << "pt=" << gammaPt
 184 //                                                      << "p=" << gammaP
 185                                                         << "theta=" << gammaTheta
 186                                                         << "coord.=" << &coord
 187                                                         << "daughterProp.=" << &daughterProp
 188                                                         << "\n";
 189                 }
 190         } // Converted MC Gamma
 191 }
 192
 193 ///________________________________________________________________________
 194 void AliAnalysisTaskMaterial::ProcessMCPhotons(){
 195         // Loop over all primary MC particle
 196         AliStack *ffMCStack = fMCEvent->Stack();
 197         for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
 198                 TParticle* particle = (TParticle *)ffMCStack->Particle(i);
 199                 if (!particle) continue;
 200                 if (particle->GetPdgCode() == 111 && particle->GetFirstDaughter() >= ffMCStack->GetNprimary()){
 201 //                      cout << "Undecayed pi0 found with mother: " << particle->GetMother(0) << endl;
 202                         for (Int_t j = 0; j < 2 ; j++){
 203                                 FillMCTree(particle->GetDaughter(j));
 204                         }
 205                 } else {
 206                         FillMCTree(i);
 207                 }
 208         }
 209 }
 210
 211 ///________________________________________________________________________
 212 void AliAnalysisTaskMaterial::ProcessPhotons(){
 213
 214         // Fill Histograms for QA and MC
 215    for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
 216       AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
 217       if (gamma ==NULL) continue;
 218       if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
 219 //              cout << "i=  " <<firstGammaIndex << " of "<< fConversionGammas->GetEntriesFast() << endl;
 220       Float_t gammaPt = gamma->GetPhotonPt();
 221 //              Float_t gammaP = gamma->GetPhotonP();
 222                 Float_t gammaTheta = gamma->GetPhotonTheta();
 223                 Float_t gammaEta = gamma->GetPhotonEta();
 224                 Float_t gammaChi2NDF = gamma->GetChi2perNDF();
 225 //              Float_t gammaX = gamma->GetConversionX();
 226 //              Float_t gammaY = gamma->GetConversionY();
 227                 Float_t gammaZ = gamma->GetConversionZ();
 228       Float_t gammaR = gamma->GetConversionRadius();
 229                 Float_t gammaPhi = gamma->GetPhotonPhi();
 230
 231                 TVectorF coord(5);
 232                 coord(0) = gamma->GetConversionX();
 233                 coord(1) = gamma->GetConversionY();
 234                 coord(2) = gamma->GetConversionZ();
 235                 coord(3) = gamma->GetConversionRadius();
 236                 coord(4) = gamma->GetPhotonPhi();
 237                 TVectorF daughterProp(4);
 238                 AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
 239       AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
 240                 daughterProp(0) = posTrack->Pt();
 241       daughterProp(1) = posTrack->Theta();
 242       daughterProp(2) = negTrack->Pt();
 243       daughterProp(3) = negTrack->Theta();
 244
 245                 UInt_t kind = 9;
 246                 if(fMCEvent){
 247 //                      cout << "generating MC stack"<< endl;
 248                         AliStack *fMCStack = fMCEvent->Stack();
 249                         if (!fMCStack) continue;
 250                         TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
 251                         TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
 252 //                      cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
 253
 254                         if(posDaughter == NULL || negDaughter == NULL){
 255                                 kind = 9; // garbage
 256 //                              cout << "one of the daughters not available" << endl;
 257                         } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
 258                                 // Not Same Mother == Combinatorial Bck
 259                                 kind = 1;
 260 //                              cout << "not the same mother" << endl;
 261                                 Int_t pdgCodePos;
 262                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 263                                 Int_t pdgCodeNeg;
 264                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 265 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 266                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
 267                                         kind = 10; //Electron Combinatorial
 268                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
 269                                         kind = 15; //direct Electron Combinatorial
 270
 271                                 if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
 272                                         kind = 11; //Pion Combinatorial
 273                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) ||
 274                                         (TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
 275                                         kind = 12; //Pion, Proton Combinatorics
 276                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) ||
 277                                         (TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
 278                                         kind = 13; //Pion, Electron Combinatorics
 279                                 if (TMath::Abs(pdgCodePos)==321 || TMath::Abs(pdgCodeNeg)==321)
 280                                         kind = 14; //Kaon combinatorics
 281
 282                         } else {
 283 //                              cout << "same mother" << endl;
 284                                 Int_t pdgCodePos;
 285                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
 286                                 Int_t pdgCodeNeg;
 287                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
 288 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
 289                                 Int_t pdgCode;
 290                                 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
 291 //                              cout << "PDG code: " << pdgCode << endl;
 292                                 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
 293                                         kind = 2; // combinatorics from hadronic decays
 294                                 else if ( !(pdgCodeNeg==pdgCodePos)){
 295                                         TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
 296                                         Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
 297                                         if(pdgCode == 111)
 298                                                 kind = 3; // pi0 Dalitz
 299                                         else if (pdgCode == 221)
 300                                                 kind = 4; // eta Dalitz
 301                                         else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
 302                                                 if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
 303                                                         kind = 0; // primary photons
 304                                                 } else if (pdgCode == 22){
 305                                                         kind = 5; //secondary photons
 306                                                 }
 307                                                 if(pdgCode == 22){
 308                                                         Float_t mcPt   = truePhotonCanditate->Pt();
 309                                                         Float_t mcR = gamma->GetNegativeMCDaughter(fMCStack)->R();
 310                                                         Float_t mcZ = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
 311                                                         Float_t mcPhi = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
 312                                                         Float_t mcEta = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
 313                                                         if (fStreamResolution){
 314                                                                 (*fStreamResolution)<<"Resolution"
 315                                                                 << "ESDpt=" << gammaPt
 316                                                                 << "ESDphi=" << gammaPhi
 317                                                                 << "ESDeta=" << gammaEta
 318                                                                 << "ESDr="<< gammaR
 319                                                                 << "ESDz="<< gammaZ
 320                                                                 << "MCpt=" << mcPt
 321                                                                 << "MCphi=" << mcPhi
 322                                                                 << "MCeta=" << mcEta
 323                                                                 << "MCr="<< mcR
 324                                                                 << "MCz="<< mcZ
 325                                                                 << "chi2ndf=" << gammaChi2NDF
 326                                                                 << "\n";
 327                                                         }
 328                                                 }
 329                                         } else  kind = 9; //garbage
 330                                 } else kind = 9; //garbage
 331                         }
 332 //                      cout << gammaPt << "\t" << gammaP<< "\t" << gammaEta<<  "\t" <<gammaChi2NDF << "\t" << gammaX<<  "\t" <<gammaY << "\t" << gammaZ<< "\t" << gammaR<<  "\t" <<  gammaPhi<< "\t" <<kind << endl;
 333
 334                         if (fStreamMaterial){
 335                                 (*fStreamMaterial)<<"ConvPointRec"
 336                                                                 << "pt=" << gammaPt
 337                                                                 << "theta=" << gammaTheta
 338                                                                 << "chi2ndf=" << gammaChi2NDF
 339                                                                 << "kind=" << kind
 340                                                                 << "coord.=" << &coord
 341                                                                 << "daughterProp.=" << &daughterProp
 342                                                                 << "\n";
 343                         }
 344                 } else {
 345                                 if (fStreamMaterial){
 346                                 (*fStreamMaterial)<<"ConvPointRec"
 347                                                                 << "pt=" << gammaPt
 348                                                                 << "theta=" << gammaTheta
 349                                                                 << "chi2ndf=" << gammaChi2NDF
 350                                                                 << "coord.=" << &coord
 351                                                                 << "daughterProp.=" << &daughterProp
 352                                                                 << "\n";
 353                         }
 354                 }
 355         }
 356 }
 357
 358 //________________________________________________________________________
 359 Int_t AliAnalysisTaskMaterial::CountESDTracks09(){
 360
 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    Int_t fNumberOfESDTracks = 0;
 369    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 370       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 371       if(!curTrack) continue;
 372       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 373    }
 374    delete EsdTrackCuts;
 375    EsdTrackCuts=0x0;
 376
 377    return fNumberOfESDTracks;
 378 }
 379
 380 Int_t AliAnalysisTaskMaterial::CountESDTracks0914(){
 381
 382    // Using standard function for setting Cuts ; We use TPCOnlyTracks for outer eta region
 383    //Bool_t selectPrimaries=kTRUE;
 384          //   EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 385    AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 386          EsdTrackCuts->SetMaxDCAToVertexXY(5);
 387          //      EsdTrackCuts->SetMaxDCAToVertexXYPtDep("sqrt(0.15^2+(0.4/pt)^2");
 388    EsdTrackCuts->SetEtaRange(0.9, 1.4);
 389    EsdTrackCuts->SetPtRange(0.15);
 390
 391    Int_t fNumberOfESDTracks = 0;
 392    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 393       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 394       if(!curTrack) continue;
 395       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 396    }
 397    EsdTrackCuts->SetEtaRange(-1.4, -0.9);
 398    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 399       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 400       if(!curTrack) continue;
 401       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 402    }
 403    delete EsdTrackCuts;
 404    EsdTrackCuts=0x0;
 405
 406    return fNumberOfESDTracks;
 407 }
 408
 409 Int_t AliAnalysisTaskMaterial::CountESDTracks14(){
 410
 411    // Using standard function for setting Cuts
 412    Bool_t selectPrimaries=kTRUE;
 413    AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 414    EsdTrackCuts->SetMaxDCAToVertexZ(2);
 415    EsdTrackCuts->SetEtaRange(-1.4, 1.4);
 416    EsdTrackCuts->SetPtRange(0.15);
 417
 418    Int_t fNumberOfESDTracks = 0;
 419    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
 420       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
 421       if(!curTrack) continue;
 422       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
 423    }
 424    delete EsdTrackCuts;
 425    EsdTrackCuts=0x0;
 426
 427    return fNumberOfESDTracks;
 428 }
 429
 430
 431 //________________________________________________________________________
 432 void AliAnalysisTaskMaterial::Terminate(Option_t *)
 433 {
 434 //    if (fStreamMaterial){
 435 //       fStreamMaterial->GetFile()->Write();
 436 //    }
 437 //    if (fStreamResolution){
 438 //       fStreamResolution->GetFile()->Write();
 439 //    }
 440 }