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