PWG1/TRD/AliTRDv0Monitor.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, 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-commercialf 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 //  Monitor V0 for TRD
  19 //
  20 //  Authors:
  21 //  Markus Heide <mheide@uni-muenster.de>
  22 //////////////////////////////////////////////////////
  23
  24 #include "TPDGCode.h"
  25 #include "TFile.h"
  26 #include "TTree.h"
  27 #include "TEventList.h"
  28 #include "TObjArray.h"
  29 #include "TH2.h"
  30 #include "TH2F.h"
  31 #include "TH1I.h"
  32 #include "TList.h"
  33 #include "TCanvas.h"
  34 #include "TPad.h"
  35 #include "TLegend.h"
  36
  37 #include "AliLog.h"
  38 #include "AliESDtrack.h"
  39 #include "AliPID.h"
  40 #include "AliTRDv0Monitor.h"
  41 #include "info/AliTRDv0Info.h"
  42
  43
  44 ClassImp(AliTRDv0Monitor)
  45
  46 //________________________________________________________________________
  47 AliTRDv0Monitor::AliTRDv0Monitor()
  48   :AliTRDrecoTask()
  49    ,fOutput(NULL)
  50    ,fhQualityReductions(NULL)
  51    ,fV0s(NULL)
  52    ,fData(NULL)
  53    ,fInfo(NULL)
  54    ,fP(-1.)
  55 {
  56   //
  57   // Default constructor
  58   //
  59   SetNameTitle("v0Monitor", "Reference V0 Monitor");
  60 }
  61
  62 //________________________________________________________________________
  63 AliTRDv0Monitor::AliTRDv0Monitor(const char *name, const char *title)
  64   :AliTRDrecoTask(name, title)
  65    ,fOutput(NULL)
  66    ,fhQualityReductions(NULL)
  67    ,fV0s(NULL)
  68    ,fData(NULL)
  69    ,fInfo(NULL)
  70    ,fP(-1.)
  71 {
  72   //
  73   // Default constructor
  74   //
  75   DefineInput(2, TObjArray::Class()); // v0 list
  76   DefineInput(3, TObjArray::Class()); // pid info list
  77   DefineOutput(1, TList::Class());
  78
  79 }
  80
  81
  82 //________________________________________________________________________
  83 AliTRDv0Monitor::~AliTRDv0Monitor()
  84 {
  85   if(fOutput) delete fOutput;
  86 }
  87 //____________________________________________________________________
  88 void AliTRDv0Monitor::MakeSummary(){
  89   TCanvas *cOut = new TCanvas("v0MonitorSummary1", "Summary 1 for task V0Monitor", 1024, 768);
  90   cOut->cd();
  91   //GetRefFigure(4);
  92   cOut->SaveAs("V0MonitorSummary.gif");
  93
  94   cOut = new TCanvas("v0MonitorSummary2","Summary 2 for task V0Monitor", 1024, 768);
  95   cOut->cd();
  96   //GetRefFigure(5);
  97   cOut->SaveAs("V0MonitorSummary2.gif");
  98 }
  99
 100 //________________________________________________________________________
 101 Bool_t AliTRDv0Monitor::GetRefFigure(Int_t /*ifig*/)
 102 {
 103
 104   AliInfo("Implementation on going ...");
 105   return kTRUE;
 106 }
 107
 108 //________________________________________________________________________
 109 void AliTRDv0Monitor::UserCreateOutputObjects()
 110 {
 111   // Create histograms
 112   // Called once
 113
 114
 115  fOutput = new TList;
 116  fOutput->SetName("V0Monitoring");
 117
 118
 119  const char *samplename[AliPID::kSPECIES] = {"electrons","muons","pions","kaons","protons"};
 120  const char *decayname[AliTRDv0Info::kNDecays] = {"gamma","K0s","Lambda","AntiLambda"};
 121  const char *detectorname[kNDets] = {"ITS","TPC","TOF"};
 122
 123  fhQualityReductions = new TH1I(Form("fhQualityReductions"),Form("Number of tracks cut out by different quality cut steps"),11,-9,2);
 124  fOutput->Add(fhQualityReductions);
 125
 126  for(Int_t ipart = 0;ipart < AliPID::kSPECIES; ipart++){
 127    fhCutReductions[ipart] = new TH1I(Form("fhCutReductions_%s",samplename[ipart]),Form("Number of tracks cut out by different cut steps for %s",samplename[ipart]),19,-17,2);
 128    fOutput->Add(fhCutReductions[ipart]);
 129    for(Int_t idetector = 0; idetector < kNDets; idetector++){
 130      fhDetPID[idetector][ipart] = new TH2F(Form("fhDetector_%s_%s",detectorname[idetector],samplename[ipart]),Form("%s Likelihood for %s vs. momentum",detectorname[idetector], samplename[ipart]),100,0.,10.,100, 0., 1.);
 131
 132      fOutput->Add(fhDetPID[idetector][ipart]);
 133    }
 134    fhComPID[ipart] = new TH2F(Form("fhComPID_%s",samplename[ipart]),Form("Combined TPC/TOF PID: Likelihood for %s",samplename[ipart]),100,0.,10.,100,0.,1.);
 135
 136      fOutput->Add(fhComPID[ipart]);
 137
 138      for(Int_t cutstep = 0; cutstep < kNCutSteps; cutstep++){
 139        fhTPCdEdx[ipart][cutstep] = new TH2F(Form("fhTPCdEdx_%s_[%d]",samplename[ipart],cutstep),Form("TPC dE/dx for %s [%d]",samplename[ipart],cutstep),100,0.,10.,300,0.,300.);
 140
 141        fOutput->Add(fhTPCdEdx[ipart][cutstep]);
 142      }
 143  }
 144
 145  for(Int_t iDecay = 0; iDecay < AliTRDv0Info::kNDecays; iDecay++){
 146    for(Int_t cutstep =0; cutstep < kNCutSteps; cutstep++){
 147      fhV0Chi2ndf[iDecay][cutstep] =  new TH2F(Form("fhV0Chi2ndf_%s_[%d]",decayname[iDecay],cutstep),Form("Chi2/NDF vs. momentum"),100,0.,10.,500, 0., 500.);
 148
 149      fOutput->Add(fhV0Chi2ndf[iDecay][cutstep]);
 150
 151      fhPsiPair[iDecay][cutstep] =  new TH2F(Form("fhV0PsiPair_%s_[%d]",decayname[iDecay],cutstep),Form("Psi_pair vs. momentum"),100,0.,10.,200, 0., 1.6);
 152
 153      fOutput->Add(fhPsiPair[iDecay][cutstep]);
 154
 155      fhPointAngle[iDecay][cutstep] =  new TH2F(Form("fhPointAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Pointing Angle vs. momentum"),100,0.,10.,500, 0., 1.6);
 156      fOutput->Add(fhPointAngle[iDecay][cutstep]);
 157
 158      fhDCA[iDecay][cutstep] =  new TH2F(Form("fhDCA_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Daughter DCA vs. momentum"),100,0.,10.,500, 0., 1.);
 159
 160      fOutput->Add(fhDCA[iDecay][cutstep]);
 161
 162      fhOpenAngle[iDecay][cutstep] =  new TH2F(Form("fhOpenAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Opening Angle vs. momentum"),100,0.,10.,500, 0., 1.6);
 163
 164      fOutput->Add(fhOpenAngle[iDecay][cutstep]);
 165
 166      fhRadius[iDecay][cutstep] =  new TH2F(Form("fhRadius_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Generation Radius vs. momentum"),100,0.,10.,500, 0., 150.);
 167
 168      fOutput->Add(fhRadius[iDecay][cutstep]);
 169
 170
 171
 172
 173    }
 174
 175    fhInvMass[iDecay] =  new TH2F(Form("fhInvMass_%s",decayname[iDecay]),Form("Invariant Mass vs. momentum"),100,0.,10.,500, 0., 2.);
 176
 177    fOutput->Add(fhInvMass[iDecay]);
 178
 179  }
 180
 181  /*TH1F *hV0mcPID[AliPID::kSPECIES][AliPID::kSPECIES];
 182    Int_t nPBins = 200;
 183
 184
 185
 186    for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++){
 187    for(Int_t iSample = 0; iSample < AliPID::kSPECIES; iSample++){
 188
 189    fhV0mcPID[iSample][iSpecies] = new TH1F(Form("fhV0mcPID_%s_is_%s",name[iSample],name[iSpecies]),Form("%s contained in %s sample",name[iSpecies],name[iSample]), nPBins, 0.2, 13.);
 190    }
 191    }*/
 192 }
 193 //________________________________________________________________________
 194 void AliTRDv0Monitor::UserExec(Option_t *)
 195 {
 196   // Main loop
 197   // Called for each event
 198
 199   if(!(fTracks = dynamic_cast<TObjArray*>(GetInputData(1)))) return;
 200   if(!(fV0s    = dynamic_cast<TObjArray*>(GetInputData(2)))) return;
 201   if(!(fInfo   = dynamic_cast<TObjArray*>(GetInputData(3)))) return;
 202
 203
 204   AliTRDtrackInfo     *track = NULL;
 205   AliTRDv0Info *v0(NULL);
 206
 207
 208
 209   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 210     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 211     for(Int_t iv(0); iv<fV0s->GetEntriesFast(); iv++){
 212       if(!(v0 = (AliTRDv0Info*)fV0s->At(iv))) continue;
 213       if(!v0->HasTrack(track)) continue;
 214       ULong_t status = track->GetStatus();
 215       if(!(status&AliESDtrack::kTRDpid)) continue;
 216
 217       fhQualityReductions->Fill(v0->fQuality);//fills integer codes for tracks cut out by track/V0 quality cuts
 218
 219       if(!(v0->fQuality == 1))continue;
 220
 221       for(Int_t part = 0; part < AliPID::kSPECIES; part++){
 222         fhCutReductions[part]->Fill(v0->GetPID(part,track));//fill in numbers of tracks eliminated by different PID cuts
 223       }
 224
 225
 226       for(Int_t idecay(0), part(-1); idecay <  AliTRDv0Info::kNDecays; idecay++){//loop over decay types considered for reference data
 227
 228         if(idecay ==  AliTRDv0Info::kLambda){ //protons and pions from Lambda
 229           part = AliPID::kProton;
 230         } else if(idecay == AliTRDv0Info::kAntiLambda) { //antiprotons and pions from Anti-Lambda
 231           part = AliPID::kProton;
 232         } else if(idecay ==   AliTRDv0Info::kK0s) {//pions from K0s
 233           part = AliPID::kPion;
 234         } else if(idecay ==  AliTRDv0Info::kGamma) {//electrons from conversions
 235           part = AliPID::kElectron;
 236         }
 237
 238         //fill histograms with track/V0 quality cuts only
 239         fhPsiPair[idecay][0]->Fill(v0->fV0Momentum,v0->fPsiPair);//Angle between daughter momentum plane and plane perpendicular to magnetic field
 240         fhInvMass[idecay]->Fill(v0->fV0Momentum,v0->fInvMass[idecay]);//Invariant mass
 241         fhPointAngle[idecay][0]->Fill(v0->fV0Momentum,v0->fPointingAngle);// = TMath::ACos(esdv0->GetV0CosineOfPointingAngle()); // Cosine of pointing angle
 242         fhOpenAngle[idecay][0]->Fill(v0->fV0Momentum,v0->fOpenAngle);// opening angle between daughters
 243         fhDCA[idecay][0]->Fill(v0->fV0Momentum,v0->fDCA);// Distance of closest approach of daughter tracks
 244         fhV0Chi2ndf[idecay][0]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);//Kalman Filter Chi2/NDF
 245         fhRadius[idecay][0]->Fill(v0->fV0Momentum,v0->fRadius);//distance of decay/conversion from primary vertex in x-y plane
 246
 247         if(v0->HasTrack(track) == -1){
 248           fhTPCdEdx[part][0]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);//TPC dE/dx for negative track
 249         }
 250         else if(v0->HasTrack(track) == 1){
 251           fhTPCdEdx[part][0]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);//TPC dE/dx for positive track
 252         }
 253         //fill histograms after invariant mass cuts
 254         if((v0->fInvMass[idecay] < v0->fUpInvMass[idecay][0])&&(v0->fInvMass[idecay]> v0->fDownInvMass[idecay])){
 255           fhV0Chi2ndf[idecay][1]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);
 256           fhPsiPair[idecay][1]->Fill(v0->fV0Momentum,v0->fPsiPair);
 257           fhPointAngle[idecay][1]->Fill(v0->fV0Momentum,v0->fPointingAngle);
 258           fhOpenAngle[idecay][1]->Fill(v0->fV0Momentum,v0->fOpenAngle);
 259           fhDCA[idecay][1]->Fill(v0->fV0Momentum,v0->fDCA);
 260           fhRadius[idecay][1]->Fill(v0->fV0Momentum,v0->fRadius);
 261           if(v0->HasTrack(track) == -1)
 262             fhTPCdEdx[part][1]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);
 263           else if(v0->HasTrack(track) == 1)
 264             fhTPCdEdx[part][1]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);
 265
 266         }
 267
 268         //fill histograms after all reference selection cuts
 269         if(v0->GetPID(part,track)==1){
 270           fhV0Chi2ndf[idecay][2]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);
 271           fhPsiPair[idecay][2]->Fill(v0->fV0Momentum,v0->fPsiPair);
 272           fhPointAngle[idecay][2]->Fill(v0->fV0Momentum,v0->fPointingAngle);
 273           fhOpenAngle[idecay][2]->Fill(v0->fV0Momentum,v0->fOpenAngle);
 274           fhDCA[idecay][2]->Fill(v0->fV0Momentum,v0->fDCA);
 275           fhRadius[idecay][2]->Fill(v0->fV0Momentum,v0->fRadius);
 276           if(v0->HasTrack(track) == -1)
 277             fhTPCdEdx[part][2]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);
 278           else if(v0->HasTrack(track) == 1)
 279             fhTPCdEdx[part][2]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);
 280
 281         }
 282       }
 283     }
 284   }
 285   PostData(1, fOutput);
 286 }
 287 //_______________________________________________________________________
 288
 289 /*Int_t AliTRDv0Monitor::GetPDG(Int_t index)
 290 {//yet to come
 291
 292
 293
 294 }*/
 295 //_______________________________________________________________________
 296