PWGDQ/dielectron/AliAnalysisTaskDielectronReadAODBranch.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   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 //  Class to retrieve the branch of the dielectron candidates stored in   //
  18 //  filtered AODs file (AliAOD.Dielectron.root). It is possible to        //
  19 //  apply tighter cuts to candidates stored in the branch and do the      //
  20 //  matching with MC truth.                                               //
  21 ////////////////////////////////////////////////////////////////////////////
  22
  23 #include "TChain.h"
  24 #include "TNtuple.h"
  25 #include "TList.h"
  26 #include "TH1F.h"
  27 #include "TH2F.h"
  28 #include "TDatabasePDG.h"
  29 #include "TF1.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliAODHandler.h"
  32 #include "AliAODEvent.h"
  33 #include "AliAODVertex.h"
  34 #include "AliAODTrack.h"
  35 #include "AliAODMCParticle.h"
  36 #include "AliAnalysisTaskDielectronReadAODBranch.h"
  37 #include "AliDielectronPair.h"
  38 #include "AliDielectronMC.h"
  39
  40 ClassImp(AliAnalysisTaskDielectronReadAODBranch)
  41         //________________________________________________________________________
  42         AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch():
  43                 AliAnalysisTaskSE(),
  44                 fOutput(0),
  45                 fNtupleJPSI(0),
  46                 fNentries(0),
  47                 fInvMass(0),
  48                 fInvMassNoCuts(0),
  49                 fpsproperSignal(0),
  50                 fpsproperSidebands(0),
  51                 fpsproperAll(0),
  52                 fpsproperUnder(0),
  53                 fpsproperUpper(0),
  54                 fLxyVsPtleg1(0),
  55                 fLxyVsPtleg2(0),
  56                 fLxyVsPt(0),
  57                 fMeeVsPt(0),
  58                 fMeeVsLxy(0),
  59                 fprimvtxZ(0),
  60                 fsecvtxZ(0),
  61                 fprimvtxX(0),
  62                 fsecvtxX(0),
  63                 fprimvtxY(0),
  64                 fsecvtxY(0),
  65                 fPt(0),
  66                 fPtLeg1(0),
  67                 fPtLeg2(0),
  68                 fdEdxP(0),
  69                 fHasMC(0),
  70                 fobj(0),
  71                 fobjMC(0),
  72                 fPtCut(1.),
  73                 fSpdFirstRequired(kFALSE),
  74                 fClsTPC(90),
  75                 fPairType(1),
  76                 fPtJpsi(1.3),
  77                 fInvMassSignalLimits(0),
  78                 fInvMassSideBandsLimits(0),
  79                 fSecondary(0)
  80 {
  81         // Default constructor
  82 }
  83
  84 //________________________________________________________________________
  85 AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch(const char *name):
  86         AliAnalysisTaskSE(name),
  87         fOutput(0),
  88         fNtupleJPSI(0),
  89         fNentries(0),
  90         fInvMass(0),
  91         fInvMassNoCuts(0),
  92         fpsproperSignal(0),
  93         fpsproperSidebands(0),
  94         fpsproperAll(0),
  95         fpsproperUnder(0),
  96         fpsproperUpper(0),
  97         fLxyVsPtleg1(0),
  98         fLxyVsPtleg2(0),
  99         fLxyVsPt(0),
 100         fMeeVsPt(0),
 101         fMeeVsLxy(0),
 102         fprimvtxZ(0),
 103         fsecvtxZ(0),
 104         fprimvtxX(0),
 105         fsecvtxX(0),
 106         fprimvtxY(0),
 107         fsecvtxY(0),
 108         fPt(0),
 109         fPtLeg1(0),
 110         fPtLeg2(0),
 111         fdEdxP(0),
 112         fHasMC(0),
 113         fobj(0),
 114         fobjMC(0),
 115         fPtCut(1.),
 116         fSpdFirstRequired(kFALSE),
 117         fClsTPC(90),
 118         fPairType(1),
 119         fPtJpsi(1.3),
 120         fInvMassSignalLimits(0),
 121         fInvMassSideBandsLimits(0),
 122         fSecondary(0)
 123 {
 124         // Default constructor
 125         DefineInput(0,TChain::Class());
 126         DefineOutput(1,TList::Class());  //My private output
 127
 128         fInvMassSignalLimits = new Double_t[2];
 129         fInvMassSignalLimits[0] = 0.;  fInvMassSignalLimits[1] = 0.;
 130         fInvMassSideBandsLimits = new Double_t[2];
 131         fInvMassSideBandsLimits[0] = 0.;  fInvMassSideBandsLimits[1] = 0.;
 132 }
 133
 134 //________________________________________________________________________
 135 AliAnalysisTaskDielectronReadAODBranch::~AliAnalysisTaskDielectronReadAODBranch()
 136 {
 137         if(fOutput)   delete fOutput;
 138         if(fobj)      delete fobj;
 139         if(fobjMC)    delete fobjMC;
 140 }
 141
 142 //________________________________________________________________________
 143 void AliAnalysisTaskDielectronReadAODBranch::Init()
 144 {
 145         // Initialization
 146         if(fDebug > 1) printf("AnalysisTaskReadAOD::Init() \n");
 147         return;
 148 }
 149
 150 //________________________________________________________________________
 151 void AliAnalysisTaskDielectronReadAODBranch::UserCreateOutputObjects()
 152 {
 153         //
 154         // Create the output container
 155         //
 156         if(fDebug > 1) printf("AnalysisTaskReadAOD::UserCreateOutputObjects() \n");
 157         fOutput = new TList();
 158         fOutput->SetOwner();
 159         //// Histogram booking
 160
 161         // invariant mass
 162         fInvMass = new TH1F("fInvMass","fInvMass",300,2.0,2.0+300*.04); // step 40MeV
 163         fInvMassNoCuts = new TH1F("fInvMass_no_cuts","fInvMass_no_cuts",125,0,125*0.04); //step 40MeV
 164         fNentries = new TH1F("numberOfevent","numbOfEvent",1,-0.5,0.5);
 165
 166         //pseudoproper
 167         fpsproperSignal = new TH1F("psproper_decay_length",Form("psproper_decay_length_distrib(AODcuts+#clustTPC>%d, pT(leg)>%fGeV, pT(J/#psi)>%fGeV/c, %f < M < %f);X [#mu m];Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]),150,-3000.,3000.);
 168         fpsproperSidebands = new TH1F("psproper_decay_length_sidebands",Form("psproper_decay_length_distrib_sidebands(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,M> %f && M < %f );  X [#mu m]; Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSideBandsLimits[1],fInvMassSideBandsLimits[0]),150,-3000.,3000.);
 169         fpsproperAll = new TH1F("psproper_decay_length_all",Form("psproper_decay_length_distrib_all(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
 170         fpsproperUnder = new TH1F("psproper_decay_length_under",Form("psproper_decay_length_distrib_under(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
 171         fpsproperUpper = new TH1F("psproper_decay_length_upper",Form("psproper_decay_length_distrib_upper(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
 172         //
 173         fLxyVsPtleg1 = new TH2F("Lxy_vs_Pt_leg1",Form("Lxy_vs_Pt_leg1_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
 174         fLxyVsPtleg2 = new TH2F("Lxy_vs_Pt_leg2",Form("Lxy_vs_Pt_leg2_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
 175         fLxyVsPt = new TH2F("Lxy_vs_Pt_jpsi",Form("Lxy_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
 176         fMeeVsPt = new TH2F("Mee_vs_Pt_jpsi",Form("Mee_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; M_{ee}",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,200, 1.99,4.1);
 177         fMeeVsLxy = new TH2F("Mee_vs_Lxy_jpsi",Form("Mee_vs_Lxy_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV),pT(J/#psi)>%f GeV/c; X; M_{ee}",fClsTPC,fPtCut,fPtJpsi), 400, -3000, 3000.,200, 1.99,4.1);
 178
 179         // QA plots
 180         fprimvtxZ = new TH1F("prim_vtx_Z",Form("prim_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f);  Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
 181         fsecvtxZ = new TH1F("sec_vtx_Z",Form("sec_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f);  Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
 182         fprimvtxX = new TH1F("prim_vtx_X",Form("prim_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  X[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
 183         fsecvtxX = new TH1F("sec_vtx_X",Form("sec_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  X[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
 184         fprimvtxY = new TH1F("prim_vtx_Y",Form("prim_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  Y[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
 185         fsecvtxY = new TH1F("sec_vtx_Y",Form("sec_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  Y[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
 186         //
 187         fPt = new TH1F("Pt(J/psi)",Form("Pt_Jpsi_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); p_{T}(J/#psi)[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
 188         fPtLeg1 = new TH1F("Pt_leg1",Form("Pt_leg1_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f);  p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
 189         fPtLeg2 = new TH1F("Pt_leg2",Form("Pt_leg2_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f);  p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
 190         //dE/dx plots
 191         fdEdxP = new TH2F("dE/dx_vs_Ptpc",Form("dE/dx_vs_Ptpc(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f)",fClsTPC,fPtCut,fPtJpsi),400,0.2,20.,200,0.,200.);
 192
 193         //J/psi n-tuple (X and M values)
 194         fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
 195
 196         // add histograms to list
 197         fOutput->Add(fNentries);
 198         fOutput->Add(fInvMass);
 199         fOutput->Add(fInvMassNoCuts);
 200         fOutput->Add(fpsproperSignal);
 201         fOutput->Add(fpsproperSidebands);
 202         fOutput->Add(fpsproperAll);
 203         fOutput->Add(fpsproperUnder);
 204         fOutput->Add(fpsproperUpper);
 205         fOutput->Add(fLxyVsPtleg1);
 206         fOutput->Add(fLxyVsPtleg2);
 207         fOutput->Add(fLxyVsPt);
 208         fOutput->Add(fMeeVsPt);
 209         fOutput->Add(fMeeVsLxy);
 210         fOutput->Add(fprimvtxZ);
 211         fOutput->Add(fsecvtxZ);
 212         fOutput->Add(fprimvtxX);
 213         fOutput->Add(fsecvtxX);
 214         fOutput->Add(fprimvtxY);
 215         fOutput->Add(fsecvtxY);
 216         fOutput->Add(fPt);
 217         fOutput->Add(fPtLeg1);
 218         fOutput->Add(fPtLeg2);
 219         fOutput->Add(fNtupleJPSI);
 220         fOutput->Add(fdEdxP);
 221         return;
 222 }
 223
 224 //________________________________________________________________________
 225 void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
 226 {
 227         // Execute analysis for current event:
 228         AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 229         if (!aod) return;
 230
 231         Double_t vtxPrim[3] = {0.,0.,0.};
 232
 233         AliAODVertex* primvtx = aod->GetPrimaryVertex();
 234         vtxPrim[0] = primvtx->GetX();
 235         vtxPrim[1] = primvtx->GetY();
 236         vtxPrim[2] = primvtx->GetZ();
 237
 238         AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 239
 240         TTree *aodTree = aodHandler->GetTree();
 241         if(!fobj) aodTree->SetBranchAddress("dielectrons",&fobj);
 242
 243         if(fHasMC){ aodTree->SetBranchAddress("mcparticles",&fobjMC);
 244         if(!fobjMC) printf("AliAnalysisTaskDielectronReadAODBranch::UserExec: MC particles branch not found!\n"); }
 245
 246         fNentries->Fill(0);
 247         aodTree->GetEvent(Entry());
 248
 249         // loop over candidates
 250         if(fobj) {
 251                 TObjArray *objArr = (TObjArray*)fobj->UncheckedAt(fPairType);
 252                 for(int j=0;j<objArr->GetEntriesFast();j++)
 253                 {
 254                         AliDielectronPair *pairObj = (AliDielectronPair*)objArr->UncheckedAt(j);
 255                         fInvMassNoCuts->Fill(pairObj->M());
 256                         Double_t vtxSec[3] = {0.,0.,0.};
 257                         fprimvtxX->Fill(vtxPrim[0]);
 258                         fprimvtxY->Fill(vtxPrim[1]);
 259                         fprimvtxZ->Fill(vtxPrim[2]);
 260                         pairObj->XvYvZv(vtxSec);
 261                         fsecvtxX->Fill(vtxSec[0]);
 262                         fsecvtxY->Fill(vtxSec[1]);
 263                         fsecvtxZ->Fill(vtxSec[2]);
 264
 265                         Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(pairObj->Px()) + (vtxSec[1]-vtxPrim[1])*(pairObj->Py()))/pairObj->Pt();
 266                         Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/pairObj->Pt();
 267
 268                         AliAODTrack *trk = (AliAODTrack*)pairObj->GetFirstDaughterP();
 269                         AliAODTrack *trk1 = (AliAODTrack*)pairObj->GetSecondDaughterP();
 270                         if(!trk || !trk1) {printf("ERROR: daughter tracks not available\n"); continue;}
 271
 272                         //check in case of MC analysis if candidate is a true J/psi->ee
 273                           if(fHasMC && ((AliDielectronMC::Instance()->IsJpsiPrimary(pairObj))!=fSecondary)) continue;
 274                            AliAODPid *pid  = trk->GetDetPid();
 275                                 AliAODPid *pid1 = trk1->GetDetPid();
 276                                 if(!pid || !pid1){
 277                                         if(fDebug>1) printf("No AliAODPid found\n");
 278                                         continue;
 279                                 }
 280                                  // ptLeg cut
 281                                 if((trk->Pt()<fPtCut) || (trk1->Pt()<fPtCut)) continue;
 282
 283                                 // pt jpsi cut
 284                                 if((pairObj->Pt())<fPtJpsi) continue;
 285
 286                               //spd first required
 287                               if(fSpdFirstRequired)
 288                               {
 289                                if((!trk->HasPointOnITSLayer(0)) || (!trk1->HasPointOnITSLayer(0))) continue;
 290                               }
 291                            //
 292                           if((trk->GetTPCNcls()) <= fClsTPC || (trk1->GetTPCNcls()) <= fClsTPC) continue;
 293                                 if(trk->Charge()==-1){
 294                                         fPtLeg1->Fill(trk->Pt());
 295                                 } else {fPtLeg2->Fill(trk->Pt());}
 296                                 if(trk1->Charge()==-1){
 297                                         fPtLeg1->Fill(trk1->Pt());
 298                                 } else {fPtLeg2->Fill(trk1->Pt());}
 299                                 //Fill dE/dx related plots (before PID cuts)
 300                                 fdEdxP->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
 301                                 fdEdxP->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
 302                                         fInvMass->Fill(pairObj->M());
 303                                         fPt->Fill(pairObj->Pt());
 304                                         fMeeVsLxy->Fill(10000*psProperDecayLength,pairObj->M());
 305                                         fMeeVsPt->Fill(pairObj->Pt(),pairObj->M());
 306                                         fpsproperAll->Fill(10000*psProperDecayLength);
 307
 308                                            //psproper in the signal region
 309                                            if((pairObj->M())<fInvMassSignalLimits[1] && (pairObj->M())>fInvMassSignalLimits[0]){
 310                                                 fpsproperSignal->Fill(10000*psProperDecayLength);
 311                                                 fLxyVsPt->Fill(pairObj->Pt(),10000*psProperDecayLength);  //jpsi
 312                                                 fLxyVsPtleg1->Fill(trk->Pt(),10000*psProperDecayLength);  //leg1 (warning: mixture of pos and neg tracks)
 313                                                 fLxyVsPtleg2->Fill(trk1->Pt(),10000*psProperDecayLength); //leg2 (warning: mixture of pos and neg tracks)
 314                                                 fNtupleJPSI->Fill(10000*psProperDecayLength,pairObj->M()); // fill the N-tuple (imput of the minimization algorithm)
 315                                                  }
 316
 317                                         // psproper in the sidebands
 318                                         if((pairObj->M())> fInvMassSideBandsLimits[1] || (pairObj->M())< fInvMassSideBandsLimits[0])
 319                                         fpsproperSidebands->Fill(10000*psProperDecayLength);
 320
 321                                         // psproper in the lower sideband
 322                                         if(pairObj->M()< fInvMassSideBandsLimits[0]) fpsproperUnder->Fill(10000*psProperDecayLength);
 323
 324                                         // psproper in the upper sideband
 325                                         if(pairObj->M()> fInvMassSideBandsLimits[1]) fpsproperUpper->Fill(10000*psProperDecayLength);
 326                         } // end loop over candidates
 327                 }// end loop over pair types
 328         // Post the data
 329         PostData(1,fOutput);
 330         return;
 331 }
 332
 333 //________________________________________________________________________
 334 void AliAnalysisTaskDielectronReadAODBranch::Terminate(Option_t */*option*/)
 335 {
 336         if(fDebug > 1) printf("AliAnalysisTaskDielectronReadAODBranch::Terminate() \n");
 337         return;
 338 }
 339