PWG2/SPECTRA/AliAnalysisTaskDedx.cxx

   1 /**************************************************************************
   2  * Author: Boris Hippolyte.                                               *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //-----------------------------------------------------------------
  15 //                     AliAnalysisTaskDedx class
  16 //            This task is for QAing the dE/dx from the ESD
  17 //              Origin: B.H. Nov2007, hippolyt@in2p3.fr
  18 //-----------------------------------------------------------------
  19
  20 //-----------------------------------------------------------------
  21 //    ---> Next upgrades
  22 //      1) extent the use of fMidPseudoRapidityFlag to primaries
  23 //      2) discuss the GetSigmaToVertex(AliESDtrack*)
  24 //      3) decide the kind of refit to be used
  25 //-----------------------------------------------------------------
  26
  27 #include "TList.h"
  28 #include "TH1F.h"
  29 #include "TH2F.h"
  30 #include "TCanvas.h"
  31 #include "TVector3.h"
  32
  33 #include "AliAnalysisTaskSE.h"
  34
  35 #include "AliESDEvent.h"
  36 #include "AliESDVertex.h"
  37 #include "AliAODEvent.h"
  38 #include "AliAODVertex.h"
  39
  40 #include "AliESDv0.h"
  41
  42 #include "AliAnalysisTaskDedx.h"
  43
  44 ClassImp(AliAnalysisTaskDedx)
  45
  46 //________________________________________________________________________
  47 AliAnalysisTaskDedx::AliAnalysisTaskDedx()
  48   : AliAnalysisTaskSE(), fAnalysisType("ESD"), fCollidingSystems(0), fListHist(),
  49     fHistPtot(0), fHistMultiplicity(0), fHistTPCDedxVsMomentum(0), fHistITSDedxVsMomentum(0),
  50     fHistMassK0(0), fHistMassLambda(0), fHistMassAntiLambda(0),
  51     fHistTPCDedxVsMomPosK0(0), fHistTPCDedxVsMomNegK0(0),
  52     fHistTPCDedxVsMomPosLambda(0), fHistTPCDedxVsMomNegLambda(0),
  53     fHistTPCDedxVsMomPosAntiLambda(0), fHistTPCDedxVsMomNegAntiLambda(0),
  54     fHistDiffInOutMomentum(0), fHistDiffPrimOutMomentum(0),
  55     fHistDiffPrimMeanMomentum(0), fHistPercPrimMeanMomentum(0), fHistPrimEta(0),
  56     fHistPercPrimMeanMomentumVsEta(0), fHistPercPrimMeanMomentumVsPrim(0),
  57
  58     fHistMultiplicityCuts(0), fHistTPCDedxVsMomentumCuts(0), fHistITSDedxVsMomentumCuts(0),
  59     fHistMassK0Cuts(0), fHistMassLambdaCuts(0), fHistMassAntiLambdaCuts(0),
  60     fHistTPCDedxVsMomPosK0Cuts(0), fHistTPCDedxVsMomNegK0Cuts(0),
  61     fHistTPCDedxVsMomPosLambdaCuts(0), fHistTPCDedxVsMomNegLambdaCuts(0),
  62     fHistTPCDedxVsMomPosAntiLambdaCuts(0), fHistTPCDedxVsMomNegAntiLambdaCuts(0),
  63     fHistDiffInOutMomentumCuts(0), fHistDiffPrimOutMomentumCuts(0),
  64     fHistDiffPrimMeanMomentumCuts(0), fHistPercPrimMeanMomentumCuts(0), fHistPrimEtaCuts(0),
  65     fHistPercPrimMeanMomentumVsEtaCuts(0), fHistPercPrimMeanMomentumVsPrimCuts(0),
  66
  67     fAllConstrainedFlag(0), fMidPseudoRapidityFlag(0),
  68
  69     fSelTrackRemoveKink(0), fSelTrackWithOnTheFlyV0(0),
  70     fSelTrackMinClustersTPC(0), fSelTrackMinClustersITS(0),
  71     fSelTrackMaxChi2PerClusterTPC(0), fSelTrackMaxChi2PerClusterITS(0),
  72     fSelTrackMaxCov11(0), fSelTrackMaxCov22(0), fSelTrackMaxCov33(0),
  73     fSelTrackMaxCov44(0), fSelTrackMaxCov55(0),
  74     fSelV0MaxDcaDaughters(0), fSelV0MinDecayLength(0)
  75 {
  76   // Dummy constructor
  77 }
  78 //________________________________________________________________________
  79 AliAnalysisTaskDedx::AliAnalysisTaskDedx(const char    *rName,
  80                                          const Bool_t   rAllConstrainedFlag,
  81                                          const Bool_t   rMidPseudoRapidityFlag,
  82                                          const Bool_t   rSelTrackRemoveKink,
  83                                          const Bool_t   rSelTrackWithOnTheFlyV0,
  84                                          const Int_t    rSelTrackMinClustersTPC,
  85                                          const Int_t    rSelTrackMinClustersITS,
  86                                          const Float_t  rSelTrackMaxChi2PerClusterTPC,
  87                                          const Float_t  rSelTrackMaxChi2PerClusterITS,
  88                                          const Double_t rSelTrackMaxCov11,
  89                                          const Double_t rSelTrackMaxCov22,
  90                                          const Double_t rSelTrackMaxCov33,
  91                                          const Double_t rSelTrackMaxCov44,
  92                                          const Double_t rSelTrackMaxCov55,
  93                                          const Double_t rSelV0MaxDcaDaughters,
  94                                          const Double_t rSelV0MinDecayLength)
  95
  96   : AliAnalysisTaskSE(rName), fAnalysisType("ESD"), fCollidingSystems(0), fListHist(),
  97     fHistPtot(0), fHistMultiplicity(0), fHistTPCDedxVsMomentum(0), fHistITSDedxVsMomentum(0),
  98     fHistMassK0(0), fHistMassLambda(0), fHistMassAntiLambda(0),
  99     fHistTPCDedxVsMomPosK0(0), fHistTPCDedxVsMomNegK0(0),
 100     fHistTPCDedxVsMomPosLambda(0), fHistTPCDedxVsMomNegLambda(0),
 101     fHistTPCDedxVsMomPosAntiLambda(0), fHistTPCDedxVsMomNegAntiLambda(0),
 102     fHistDiffInOutMomentum(0), fHistDiffPrimOutMomentum(0),
 103     fHistDiffPrimMeanMomentum(0), fHistPercPrimMeanMomentum(0), fHistPrimEta(0),
 104     fHistPercPrimMeanMomentumVsEta(0), fHistPercPrimMeanMomentumVsPrim(0),
 105
 106     fHistMultiplicityCuts(0), fHistTPCDedxVsMomentumCuts(0), fHistITSDedxVsMomentumCuts(0),
 107     fHistMassK0Cuts(0), fHistMassLambdaCuts(0), fHistMassAntiLambdaCuts(0),
 108     fHistTPCDedxVsMomPosK0Cuts(0), fHistTPCDedxVsMomNegK0Cuts(0),
 109     fHistTPCDedxVsMomPosLambdaCuts(0), fHistTPCDedxVsMomNegLambdaCuts(0),
 110     fHistTPCDedxVsMomPosAntiLambdaCuts(0), fHistTPCDedxVsMomNegAntiLambdaCuts(0),
 111     fHistDiffInOutMomentumCuts(0), fHistDiffPrimOutMomentumCuts(0),
 112     fHistDiffPrimMeanMomentumCuts(0), fHistPercPrimMeanMomentumCuts(0), fHistPrimEtaCuts(0),
 113     fHistPercPrimMeanMomentumVsEtaCuts(0), fHistPercPrimMeanMomentumVsPrimCuts(0),
 114
 115     fAllConstrainedFlag(rAllConstrainedFlag), fMidPseudoRapidityFlag(rMidPseudoRapidityFlag),
 116
 117     fSelTrackRemoveKink(rSelTrackRemoveKink), fSelTrackWithOnTheFlyV0(rSelTrackWithOnTheFlyV0),
 118     fSelTrackMinClustersTPC(rSelTrackMinClustersTPC), fSelTrackMinClustersITS(rSelTrackMinClustersITS),
 119     fSelTrackMaxChi2PerClusterTPC(rSelTrackMaxChi2PerClusterTPC), fSelTrackMaxChi2PerClusterITS(rSelTrackMaxChi2PerClusterITS),
 120     fSelTrackMaxCov11(rSelTrackMaxCov11), fSelTrackMaxCov22(rSelTrackMaxCov22), fSelTrackMaxCov33(rSelTrackMaxCov33),
 121     fSelTrackMaxCov44(rSelTrackMaxCov44), fSelTrackMaxCov55(rSelTrackMaxCov55),
 122     fSelV0MaxDcaDaughters(rSelV0MaxDcaDaughters), fSelV0MinDecayLength(rSelV0MinDecayLength)
 123 {
 124   // Constructor
 125   // Define output slots only here
 126   // Output slot #1 writes into a TList container
 127   DefineOutput(1, TList::Class());
 128 }
 129 //________________________________________________________________________
 130 void AliAnalysisTaskDedx::UserCreateOutputObjects()
 131 {
 132   // Create histograms
 133   // Called once
 134
 135   fListHist = new TList();
 136   if (!fHistPtot) {
 137     fHistPtot = new TH1F("fHistPtot", "P_{tot} distribution;P_{tot} (GeV/c);dN/dP_{tot} (c/GeV)", 15, 0.1, 3.1);
 138     fHistPtot->SetMarkerStyle(kFullCircle);
 139     fListHist->Add(fHistPtot);
 140   }
 141   if (!fHistMultiplicity) {
 142     if (fCollidingSystems)
 143       fHistMultiplicity = new TH1F("fHistMultiplicity", "Multiplicity distribution;Number of tracks;Events", 200, 0, 40000);
 144     else
 145       fHistMultiplicity = new TH1F("fHistMultiplicity", "Multiplicity distribution;Number of tracks;Events", 250, 0, 250);
 146     fHistMultiplicity->SetMarkerStyle(kFullCircle);
 147     fListHist->Add(fHistMultiplicity);
 148   }
 149   if (!fHistTPCDedxVsMomentum) {
 150     fHistTPCDedxVsMomentum = new TH2F("h2TPCDedxVsMomentum","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 151     fListHist->Add(fHistTPCDedxVsMomentum);
 152   }
 153   if (!fHistITSDedxVsMomentum) {
 154     fHistITSDedxVsMomentum = new TH2F("h2ITSDedxVsMomentum","Bethe-Bloch Distribution for ITS;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 155     fListHist->Add(fHistITSDedxVsMomentum);
 156   }
 157   if (!fHistMassK0) {
 158     fHistMassK0 = new TH1F("h1MassK0","K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts",100,0.4,0.6);
 159     fListHist->Add(fHistMassK0);
 160   }
 161   if (!fHistMassLambda) {
 162     fHistMassLambda = new TH1F("h1MassLambda","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts",75,1.05,1.2);
 163     fListHist->Add(fHistMassLambda);
 164   }
 165   if (!fHistMassAntiLambda) {
 166     fHistMassAntiLambda = new TH1F("h1MassAntiLambda","#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts",75,1.05,1.2);
 167     fListHist->Add(fHistMassAntiLambda);
 168   }
 169   if (!fHistTPCDedxVsMomPosK0) {
 170     fHistTPCDedxVsMomPosK0 = new TH2F("h2TPCDedxVsMomPosK0","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 171     fListHist->Add(fHistTPCDedxVsMomPosK0);
 172   }
 173   if (!fHistTPCDedxVsMomNegK0) {
 174     fHistTPCDedxVsMomNegK0 = new TH2F("h2TPCDedxVsMomNegK0","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 175     fListHist->Add(fHistTPCDedxVsMomNegK0);
 176   }
 177   if (!fHistTPCDedxVsMomPosLambda) {
 178     fHistTPCDedxVsMomPosLambda = new TH2F("h2TPCDedxVsMomPosLambda","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 179     fListHist->Add(fHistTPCDedxVsMomPosLambda);
 180   }
 181   if (!fHistTPCDedxVsMomNegLambda) {
 182     fHistTPCDedxVsMomNegLambda = new TH2F("h2TPCDedxVsMomNegLambda","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 183     fListHist->Add(fHistTPCDedxVsMomNegLambda);
 184   }
 185   if (!fHistTPCDedxVsMomPosAntiLambda) {
 186     fHistTPCDedxVsMomPosAntiLambda = new TH2F("h2TPCDedxVsMomPosAntiLambda","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 187     fListHist->Add(fHistTPCDedxVsMomPosAntiLambda);
 188   }
 189   if (!fHistTPCDedxVsMomNegAntiLambda) {
 190     fHistTPCDedxVsMomNegAntiLambda = new TH2F("h2TPCDedxVsMomNegAntiLambda","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 191     fListHist->Add(fHistTPCDedxVsMomNegAntiLambda);
 192   }
 193   if (!fHistDiffInOutMomentum) {
 194     fHistDiffInOutMomentum = new TH1F("h1DiffInOutMomentum","Momentum Difference In-Out;Momentum (GeV/c);Counts",200,-0.1,0.1);
 195     fListHist->Add(fHistDiffInOutMomentum);
 196   }
 197   if (!fHistDiffPrimOutMomentum) {
 198     fHistDiffPrimOutMomentum = new TH1F("h1DiffPrimOutMomentum","Momentum Difference Prim-Out;Momentum (GeV/c);Counts",200,-0.1,0.1);
 199     fListHist->Add(fHistDiffPrimOutMomentum);
 200   }
 201   if (!fHistDiffPrimMeanMomentum) {
 202     fHistDiffPrimMeanMomentum = new TH1F("h1DiffPrimMeanMomentum","Momentum Difference Prim-Mean;Momentum (GeV/c);Counts",200,-0.1,0.1);
 203     fListHist->Add(fHistDiffPrimMeanMomentum);
 204   }
 205   if (!fHistPercPrimMeanMomentum) {
 206     fHistPercPrimMeanMomentum = new TH1F("h1PercPrimMeanMomentum","Momentum Percentage (Prim-Mean)/Prim;(Primary-Mean)/Primary (%);Counts",200,-10,10);
 207     fListHist->Add(fHistPercPrimMeanMomentum);
 208   }
 209   if (!fHistPrimEta) {
 210     fHistPrimEta = new TH1F("h1PrimEta","Pseudorapidity Distribution (Primaries);#eta;Counts",300,-1.5,1.5);
 211     fListHist->Add(fHistPrimEta);
 212   }
 213   if (!fHistPercPrimMeanMomentumVsEta) {
 214     fHistPercPrimMeanMomentumVsEta = new TH2F("h2PercPrimMeanMomentumVsEta","Momentum Percentage (Prim-Mean)/Prim vs #eta;#eta;(Primary-Mean)/Primary (%)",100,-1.5,1.5,100,-10,10);
 215     fListHist->Add(fHistPercPrimMeanMomentumVsEta);
 216   }
 217   if (!fHistPercPrimMeanMomentumVsPrim) {
 218     fHistPercPrimMeanMomentumVsPrim = new TH2F("h2PercPrimMeanMomentumVsPrim","Momentum Percentage (Prim-Mean)/Prim vs Prim;Momentum (GeV/c);(Primary-Mean)/Primary (%)",1500,0,15,100,-10,10);
 219     fListHist->Add(fHistPercPrimMeanMomentumVsPrim);
 220   }
 221   // Histograms after selections
 222   if (!fHistMultiplicityCuts) {
 223     if (fCollidingSystems)
 224       fHistMultiplicityCuts = new TH1F("fHistMultiplicityCuts", "Multiplicity distribution;Number of tracks;Events", 200, 0, 40000);
 225     else
 226       fHistMultiplicityCuts = new TH1F("h1HistMultiplicityCuts","Number of Good tracks;Number of tracks;Event", 250, 0, 250);
 227     fListHist->Add(fHistMultiplicityCuts);
 228   }
 229   if (!fHistTPCDedxVsMomentumCuts) {
 230     fHistTPCDedxVsMomentumCuts = new TH2F("h2TPCDedxVsMomentumCuts","Bethe-Bloch Distribution for TPC w/cuts;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 231     fListHist->Add(fHistTPCDedxVsMomentumCuts);
 232   }
 233   if (!fHistITSDedxVsMomentumCuts) {
 234     fHistITSDedxVsMomentumCuts = new TH2F("h2ITSDedxVsMomentumCuts","Bethe-Bloch Distribution for ITS w/cuts;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 235     fListHist->Add(fHistITSDedxVsMomentumCuts);
 236   }
 237   if (!fHistMassK0Cuts) {
 238     fHistMassK0Cuts = new TH1F("h1MassK0Cuts","K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts",100,0.4,0.6);
 239     fListHist->Add(fHistMassK0Cuts);
 240   }
 241   if (!fHistMassLambdaCuts) {
 242     fHistMassLambdaCuts = new TH1F("h1MassLambdaCuts","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts",75,1.05,1.2);
 243     fListHist->Add(fHistMassLambdaCuts);
 244   }
 245   if (!fHistMassAntiLambdaCuts) {
 246     fHistMassAntiLambdaCuts = new TH1F("h1MassAntiLambdaCuts","#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts",75,1.05,1.2);
 247     fListHist->Add(fHistMassAntiLambdaCuts);
 248   }
 249   if (!fHistTPCDedxVsMomPosK0Cuts) {
 250     fHistTPCDedxVsMomPosK0Cuts = new TH2F("h2TPCDedxVsMomPosK0Cuts","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 251     fListHist->Add(fHistTPCDedxVsMomPosK0Cuts);
 252   }
 253   if (!fHistTPCDedxVsMomNegK0Cuts) {
 254     fHistTPCDedxVsMomNegK0Cuts = new TH2F("h2TPCDedxVsMomNegK0Cuts","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 255     fListHist->Add(fHistTPCDedxVsMomNegK0Cuts);
 256   }
 257   if (!fHistTPCDedxVsMomPosLambdaCuts) {
 258     fHistTPCDedxVsMomPosLambdaCuts = new TH2F("h2TPCDedxVsMomPosLambdaCuts","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 259     fListHist->Add(fHistTPCDedxVsMomPosLambdaCuts);
 260   }
 261   if (!fHistTPCDedxVsMomNegLambdaCuts) {
 262     fHistTPCDedxVsMomNegLambdaCuts = new TH2F("h2TPCDedxVsMomNegLambdaCuts","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 263     fListHist->Add(fHistTPCDedxVsMomNegLambdaCuts);
 264   }
 265   if (!fHistTPCDedxVsMomPosAntiLambdaCuts) {
 266     fHistTPCDedxVsMomPosAntiLambdaCuts = new TH2F("h2TPCDedxVsMomPosAntiLambdaCuts","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 267     fListHist->Add(fHistTPCDedxVsMomPosAntiLambdaCuts);
 268   }
 269   if (!fHistTPCDedxVsMomNegAntiLambdaCuts) {
 270     fHistTPCDedxVsMomNegAntiLambdaCuts = new TH2F("h2TPCDedxVsMomNegAntiLambdaCuts","Bethe-Bloch Distribution for TPC;Momentum (GeV/c);dE/dx (A.U.)",1500,0,15,100,0,100);
 271     fListHist->Add(fHistTPCDedxVsMomNegAntiLambdaCuts);
 272   }
 273   if (!fHistDiffInOutMomentumCuts) {
 274     fHistDiffInOutMomentumCuts = new TH1F("h1DiffInOutMomentumCuts","Momentum Difference In-Out;Momentum (GeV/c);Counts",200,-0.1,0.1);
 275     fListHist->Add(fHistDiffInOutMomentumCuts);
 276   }
 277   if (!fHistDiffPrimOutMomentumCuts) {
 278     fHistDiffPrimOutMomentumCuts = new TH1F("h1DiffPrimOutMomentumCuts","Momentum Difference Prim-Out;Momentum (GeV/c);Counts",200,-0.1,0.1);
 279     fListHist->Add(fHistDiffPrimOutMomentumCuts);
 280   }
 281   if (!fHistDiffPrimMeanMomentumCuts) {
 282     fHistDiffPrimMeanMomentumCuts = new TH1F("h1DiffPrimMeanMomentumCuts","Momentum Difference Prim-Mean;Momentum (GeV/c);Counts",200,-0.1,0.1);
 283     fListHist->Add(fHistDiffPrimMeanMomentumCuts);
 284   }
 285   if (!fHistPercPrimMeanMomentumCuts) {
 286     fHistPercPrimMeanMomentumCuts = new TH1F("h1PercPrimMeanMomentumCuts","Momentum Percentage (Prim-Mean)/Prim;(Primary-Mean)/Primary (%);Counts",200,-10,10);
 287     fListHist->Add(fHistPercPrimMeanMomentumCuts);
 288   }
 289   if (!fHistPrimEtaCuts) {
 290     fHistPrimEtaCuts = new TH1F("h1PrimEtaCuts","Pseudorapidity Distribution (Primaries);#eta;Counts",300,-1.5,1.5);
 291     fListHist->Add(fHistPrimEtaCuts);
 292   }
 293   if (!fHistPercPrimMeanMomentumVsEtaCuts) {
 294     fHistPercPrimMeanMomentumVsEtaCuts = new TH2F("h2PercPrimMeanMomentumVsEtaCuts","Momentum Percentage (Prim-Mean)/Prim vs #eta;#eta;(Primary-Mean)/Primary (%)",100,-1.5,1.5,100,-10,10);
 295     fListHist->Add(fHistPercPrimMeanMomentumVsEtaCuts);
 296   }
 297   if (!fHistPercPrimMeanMomentumVsPrimCuts) {
 298     fHistPercPrimMeanMomentumVsPrimCuts = new TH2F("h2PercPrimMeanMomentumVsPrimCuts","Momentum Percentage (Prim-Mean)/Prim vs Prim;Momentum (GeV/c);(Primary-Mean)/Primary (%)",1500,0,15,100,-10,10);
 299     fListHist->Add(fHistPercPrimMeanMomentumVsPrimCuts);
 300   }
 301 }
 302
 303 //________________________________________________________________________
 304 void AliAnalysisTaskDedx::UserExec(Option_t *)
 305 {
 306   // Main loop
 307   // Called for each event
 308   AliVEvent* lEvent = InputEvent();
 309   if (!lEvent) {
 310     Printf("ERROR: Event not available");
 311     return;
 312   }
 313   Int_t    ntracks        = lEvent->GetNumberOfTracks();
 314   Printf("dE/dx analysis task: There are %d tracks in this event", ntracks);
 315
 316   Int_t    nAllTracks        = 0, nGoodTracks       = 0;
 317   Double_t lTPCDedx          = 0, lITSDedx          = 0;
 318
 319   if(fAnalysisType == "ESD") {
 320
 321   // Track loop
 322   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 323     AliESDtrack* track = ((AliESDEvent*)lEvent)->GetTrack(iTracks);
 324     if (!track) {
 325       Printf("ERROR: Could not retrieve track %d", iTracks);
 326       continue;
 327     }
 328     fHistPtot->Fill(track->P()); // Fill a ptot spectrum before any selection
 329
 330     Double_t lPrimvtxMomemtum[3];
 331     Bool_t lIsTrackConstrained = track->GetConstrainedPxPyPz(lPrimvtxMomemtum);
 332     if (fAllConstrainedFlag && (!lIsTrackConstrained)) continue; // Constrained or not to Prim Vertex.
 333
 334     const AliExternalTrackParam *lInnerETP = track->GetInnerParam();
 335     const AliExternalTrackParam *lOuterETP = track->GetOuterParam();
 336     if ((!lInnerETP)||(!lOuterETP)) continue;  // No inner nor outer params for this track.
 337
 338     Double_t lInnerMomentum = 0, lOuterMomentum = 0, lMeanMomentum = 0;
 339     if (lInnerETP) lInnerMomentum = lInnerETP->GetP(); // The "inner" momentum.
 340     if (lOuterETP) lOuterMomentum = lOuterETP->GetP(); // The "outer" momentum.
 341     if (lInnerETP&&lOuterETP) lMeanMomentum = (lInnerMomentum+lOuterMomentum)/2.; // The "mean" momentum
 342
 343     lTPCDedx = track->GetTPCsignal();
 344     lITSDedx = track->GetITSsignal();
 345     // This needs to be discussed and implemented
 346     // nSigmaToVertex = GetSigmaToVertex(track);
 347
 348     Float_t lPrimPseudoRap = 0, lPrimMomentum = 0;
 349     if (lIsTrackConstrained){
 350       TVector3 pPrimvtx(lPrimvtxMomemtum);
 351       lPrimPseudoRap = pPrimvtx.Eta();
 352       lPrimMomentum = pPrimvtx.Mag();
 353     }
 354
 355     if (
 356         ((!fAllConstrainedFlag)||(TMath::Abs(lPrimPseudoRap) < 0.5)) &&
 357         (lTPCDedx)
 358         )
 359       // fAllConstrained = 0 -> !fAllConstrained = 1 so ok anyway
 360       // fAllConstrained = 1 -> !fAllConstrained = 0 so it depends on eta condition
 361       {
 362         nAllTracks++;
 363         fHistTPCDedxVsMomentum->Fill(lInnerMomentum,lTPCDedx);
 364         fHistDiffInOutMomentum->Fill(lInnerMomentum-lOuterMomentum);
 365         if (lIsTrackConstrained){
 366           fHistITSDedxVsMomentum->Fill(lPrimMomentum,lITSDedx);
 367           fHistDiffPrimOutMomentum->Fill(lPrimMomentum-lOuterMomentum);
 368           fHistDiffPrimMeanMomentum->Fill(lPrimMomentum-lMeanMomentum);
 369           if (lPrimMomentum){
 370             fHistPercPrimMeanMomentum->Fill(100.*(lPrimMomentum-lMeanMomentum)/lPrimMomentum);
 371             fHistPercPrimMeanMomentumVsEta->Fill(lPrimPseudoRap,100.*(lPrimMomentum-lMeanMomentum)/lPrimMomentum);
 372             fHistPercPrimMeanMomentumVsPrim->Fill(lPrimMomentum,100.*(lPrimMomentum-lMeanMomentum)/lPrimMomentum);
 373           }
 374           fHistPrimEta->Fill(lPrimPseudoRap);
 375         }
 376         if(
 377            // (!fAllConstrainedFlag || (nSigmaToVertex < 3)) &&
 378            (IsAccepted(track))
 379            ){
 380           fHistTPCDedxVsMomentumCuts->Fill(lInnerMomentum,lTPCDedx);
 381           fHistDiffInOutMomentumCuts->Fill(lInnerMomentum-lOuterMomentum);
 382           if (lIsTrackConstrained){
 383             fHistITSDedxVsMomentumCuts->Fill(lPrimMomentum,lITSDedx);
 384             fHistDiffPrimOutMomentumCuts->Fill(lPrimMomentum-lOuterMomentum);
 385             fHistDiffPrimMeanMomentumCuts->Fill(lPrimMomentum-lMeanMomentum);
 386             if (lPrimMomentum){
 387               fHistPercPrimMeanMomentumCuts->Fill(100.*(lPrimMomentum-lMeanMomentum)/lPrimMomentum);
 388               fHistPercPrimMeanMomentumVsEtaCuts->Fill(lPrimPseudoRap,100.*(lPrimMomentum-lMeanMomentum)/lPrimMomentum);
 389               fHistPercPrimMeanMomentumVsPrimCuts->Fill(lPrimMomentum,100.*(lPrimMomentum-lMeanMomentum)/lPrimMomentum);
 390             }
 391             fHistPrimEtaCuts->Fill(lPrimPseudoRap);
 392           }
 393           nGoodTracks++;
 394         }
 395       }
 396   } //track loop
 397   fHistMultiplicity->Fill(ntracks);
 398   fHistMultiplicityCuts->Fill(nGoodTracks);
 399
 400   Int_t nv0s = 0;
 401   nv0s = lEvent->GetNumberOfV0s();
 402
 403   Int_t    lIndexTrackPos       = 0,  lIndexTrackNeg       = 0;
 404   Double_t lTPCDedxPos          = 0,  lTPCDedxNeg          = 0;
 405   Float_t lInvMassK0 = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
 406   for (Int_t iV0 = 0; iV0 < nv0s; iV0++)
 407     {// This is the V0 loop
 408       AliESDv0 *v0 = ((AliESDEvent*)lEvent)->GetV0(iV0);
 409       if (!v0) continue;
 410       if (v0->GetOnFlyStatus() != fSelTrackWithOnTheFlyV0) continue;
 411
 412       const AliESDVertex *lPrimaryVertex = ((AliESDEvent*)lEvent)->GetPrimaryVertex();
 413       Double_t tPositionPrimaryVertex[3]; lPrimaryVertex->GetXYZ(tPositionPrimaryVertex);
 414       Double_t v0PositionX = 0, v0PositionY = 0, v0PositionZ = 0; v0->GetXYZ(v0PositionX,v0PositionY,v0PositionZ);
 415       Double_t v0DecayLength = TMath::Sqrt(
 416                                            (v0PositionX-tPositionPrimaryVertex[0])*(v0PositionX-tPositionPrimaryVertex[0])+
 417                                            (v0PositionY-tPositionPrimaryVertex[1])*(v0PositionY-tPositionPrimaryVertex[1])+
 418                                            (v0PositionZ-tPositionPrimaryVertex[2])*(v0PositionZ-tPositionPrimaryVertex[2])
 419                                            );
 420       if (v0->GetDcaV0Daughters() > fSelV0MaxDcaDaughters) continue;
 421       if (v0DecayLength < fSelV0MinDecayLength) continue;
 422
 423       // Getting invariant mass infos directly from ESD
 424       v0->ChangeMassHypothesis(310);
 425       lInvMassK0 = v0->GetEffMass();
 426       v0->ChangeMassHypothesis(3122);
 427       lInvMassLambda = v0->GetEffMass();
 428       v0->ChangeMassHypothesis(-3122);
 429       lInvMassAntiLambda = v0->GetEffMass();
 430
 431       // Filling invariant mass histos for all candidates
 432       fHistMassK0->Fill(lInvMassK0);
 433       fHistMassLambda->Fill(lInvMassLambda);
 434       fHistMassAntiLambda->Fill(lInvMassAntiLambda);
 435
 436       // Accessing the daughter track infos
 437       lIndexTrackPos = TMath::Abs(v0->GetPindex());
 438       lIndexTrackNeg = TMath::Abs(v0->GetNindex());
 439       AliESDtrack* trackPos =((AliESDEvent*)lEvent)->GetTrack(lIndexTrackPos);
 440       AliESDtrack* trackNeg =((AliESDEvent*)lEvent)->GetTrack(lIndexTrackNeg);
 441
 442       Double_t lPrimvtxMomPos[3];
 443       Bool_t lIsPosConstrained = trackPos->GetConstrainedPxPyPz(lPrimvtxMomPos);
 444       if (fAllConstrainedFlag && (!lIsPosConstrained)) continue; // Constrained or not to Prim Vertex.
 445
 446       Double_t lPrimvtxMomNeg[3];
 447       Bool_t lIsNegConstrained = trackNeg->GetConstrainedPxPyPz(lPrimvtxMomNeg);
 448       if (fAllConstrainedFlag && (!lIsNegConstrained)) continue; // Constrained or not to Prim Vertex.
 449
 450       const AliExternalTrackParam *lInnerPosETP = trackPos->GetInnerParam();
 451       const AliExternalTrackParam *lInnerNegETP = trackNeg->GetInnerParam();
 452       if ((!lInnerPosETP)||(!lInnerNegETP)) continue; // No inner params for at least one track.
 453
 454       Double_t lInnerMomPos = 0, lInnerMomNeg = 0;
 455       if (lInnerPosETP) lInnerMomPos = lInnerPosETP->GetP();
 456       if (lInnerNegETP) lInnerMomNeg = lInnerNegETP->GetP();
 457
 458       lTPCDedxPos = trackPos->GetTPCsignal();
 459       lTPCDedxNeg = trackNeg->GetTPCsignal();
 460
 461       Float_t lPrimPseudoRapPos = 0, lPrimPseudoRapNeg = 0;
 462       if (lIsPosConstrained && lIsNegConstrained){
 463         TVector3 pPrimPos(lPrimvtxMomPos);
 464         TVector3 pPrimNeg(lPrimvtxMomNeg);
 465         lPrimPseudoRapPos = pPrimPos.Eta();
 466         lPrimPseudoRapNeg = pPrimNeg.Eta();
 467       }
 468
 469       if (TMath::Abs(lInvMassK0-0.497)<0.01) {
 470         fHistMassK0Cuts->Fill(lInvMassK0);
 471         fHistTPCDedxVsMomPosK0->Fill(lInnerMomPos,lTPCDedxPos);
 472         fHistTPCDedxVsMomNegK0->Fill(lInnerMomNeg,lTPCDedxNeg);
 473         if (
 474             (!fMidPseudoRapidityFlag) || ( TMath::Abs(lPrimPseudoRapPos) < 0.5) &&
 475             (IsAccepted(trackPos))
 476             )
 477           fHistTPCDedxVsMomPosK0Cuts->Fill(lInnerMomPos,lTPCDedxPos);
 478         if (
 479             (!fMidPseudoRapidityFlag) || ( TMath::Abs(lPrimPseudoRapNeg) < 0.5) &&
 480             (IsAccepted(trackNeg))
 481             )
 482           fHistTPCDedxVsMomNegK0Cuts->Fill(lInnerMomNeg,lTPCDedxNeg);
 483       }
 484       if (TMath::Abs(lInvMassLambda-1.115)<0.01) {
 485         fHistMassLambdaCuts->Fill(lInvMassLambda);
 486         fHistTPCDedxVsMomPosLambda->Fill(lInnerMomPos,lTPCDedxPos);
 487         fHistTPCDedxVsMomNegLambda->Fill(lInnerMomNeg,lTPCDedxNeg);
 488         if (
 489             (!fMidPseudoRapidityFlag) || ( TMath::Abs(lPrimPseudoRapPos) < 0.5) &&
 490             (IsAccepted(trackPos))
 491             )
 492           fHistTPCDedxVsMomPosLambdaCuts->Fill(lInnerMomPos,lTPCDedxPos);
 493         if (
 494             (!fMidPseudoRapidityFlag) || ( TMath::Abs(lPrimPseudoRapNeg) < 0.5) &&
 495             (IsAccepted(trackNeg))
 496             )
 497           fHistTPCDedxVsMomNegLambdaCuts->Fill(lInnerMomNeg,lTPCDedxNeg);
 498       }
 499       if (TMath::Abs(lInvMassAntiLambda-1.115)<0.01) {
 500         fHistMassAntiLambdaCuts->Fill(lInvMassAntiLambda);
 501         fHistTPCDedxVsMomPosAntiLambda->Fill(lInnerMomPos,lTPCDedxPos);
 502         fHistTPCDedxVsMomNegAntiLambda->Fill(lInnerMomNeg,lTPCDedxNeg);
 503         if (
 504             (!fMidPseudoRapidityFlag) || ( TMath::Abs(lPrimPseudoRapPos) < 0.5) &&
 505             (IsAccepted(trackPos))
 506             )
 507           fHistTPCDedxVsMomPosAntiLambdaCuts->Fill(lInnerMomPos,lTPCDedxPos);
 508         if (
 509             (!fMidPseudoRapidityFlag) || ( TMath::Abs(lPrimPseudoRapNeg) < 0.5) &&
 510             (IsAccepted(trackNeg))
 511             )
 512           fHistTPCDedxVsMomNegAntiLambdaCuts->Fill(lInnerMomNeg,lTPCDedxNeg);
 513       }
 514     }//V0 loop
 515
 516   }
 517   // Post output data.
 518   PostData(1, fListHist);
 519 }
 520
 521 //________________________________________________________________________
 522 void AliAnalysisTaskDedx::Terminate(Option_t *)
 523 {
 524   // Draw result to the screen
 525   // Called once at the end of the query
 526
 527   Printf("Reminder for options and selections: fAllConstrainedFlag = %d fMidPseudoRapidityFlag = %d fSelTrackRemoveKink = %d fSelTrackWithOnTheFlyV0 = %d fSelTrackMinClustersTPC = %d fSelTrackMinClustersITS = %d fSelTrackMaxChi2PerClusterTPC = %.2f fSelTrackMaxChi2PerClusterITS = %.2f fSelTrackMaxCov11 = %.2f fSelTrackMaxCov22 = %.2f fSelTrackMaxCov33 = %.2f fSelTrackMaxCov44 = %.2f fSelTrackMaxCov55 = %.2f fSelV0MaxDcaDaughters = %.2f fSelV0MinDecayLength= %.2f",
 528          fAllConstrainedFlag,
 529          fMidPseudoRapidityFlag,fSelTrackRemoveKink,fSelTrackWithOnTheFlyV0,
 530          fSelTrackMinClustersTPC,fSelTrackMinClustersITS,
 531          fSelTrackMaxChi2PerClusterTPC,fSelTrackMaxChi2PerClusterITS,
 532          fSelTrackMaxCov11,fSelTrackMaxCov22,fSelTrackMaxCov33,fSelTrackMaxCov44,fSelTrackMaxCov55,
 533          fSelV0MaxDcaDaughters,fSelV0MinDecayLength);
 534
 535   fHistPtot = dynamic_cast<TH1F*> (((TList*)GetOutputData(1))->FindObject("fHistPtot"));
 536   if (!fHistPtot) {
 537     Printf("ERROR: fHistPtot not available");
 538     return;
 539   }
 540
 541   TCanvas *cTaskDedx = new TCanvas("AliAnalysisTaskESDDedx","Ptot",10,10,510,510);
 542   if(fHistPtot->GetMaximum() > 0.) cTaskDedx->cd(1)->SetLogy();
 543   fHistPtot->DrawCopy("E");
 544 }
 545 //____________________________________________________________________//
 546 Bool_t AliAnalysisTaskDedx::IsAccepted(AliESDtrack* track) {
 547   // Checks if the track is excluded from the cuts
 548
 549   Int_t   fIdxInt[200];
 550   Int_t   nClustersITS = track->GetITSclusters(fIdxInt);
 551   Int_t   nClustersTPC = track->GetTPCclusters(fIdxInt);
 552
 553   Float_t chi2PerClusterITS = -1;
 554   Float_t chi2PerClusterTPC = -1;
 555   if (nClustersTPC!=0)
 556     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 557
 558   Double_t extCov[15];
 559   track->GetExternalCovariance(extCov);
 560
 561   if (fSelTrackRemoveKink)
 562     if (track->GetKinkIndex(0)>0) return kFALSE;
 563   if (fSelTrackMinClustersTPC)
 564     if (nClustersTPC < fSelTrackMinClustersTPC) return kFALSE;
 565   if (fSelTrackMinClustersITS)
 566     if (nClustersITS < fSelTrackMinClustersITS) return kFALSE;
 567   if (fSelTrackMaxChi2PerClusterTPC)
 568     if (chi2PerClusterTPC > fSelTrackMaxChi2PerClusterTPC) return kFALSE;
 569   if (fSelTrackMaxChi2PerClusterITS)
 570     if (chi2PerClusterITS > fSelTrackMaxChi2PerClusterITS) return kFALSE;
 571
 572   if (fSelTrackMaxCov11)
 573     if (extCov[0]  > fSelTrackMaxCov11) return kFALSE;
 574   if (fSelTrackMaxCov22)
 575     if (extCov[2]  > fSelTrackMaxCov22) return kFALSE;
 576   if (fSelTrackMaxCov33)
 577     if (extCov[5]  > fSelTrackMaxCov33) return kFALSE;
 578   if (fSelTrackMaxCov44)
 579     if (extCov[9]  > fSelTrackMaxCov44) return kFALSE;
 580   if(fSelTrackMaxCov55)
 581     if (extCov[14] > fSelTrackMaxCov55) return kFALSE;
 582   /*
 583   if(fMaxSigmaToVertexFlag)
 584     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
 585   if(fITSRefitFlag)
 586     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
 587   if(fTPCRefitFlag)
 588     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
 589   if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
 590   */
 591   return kTRUE;
 592 }