TRD/qaAnalysis/AliTRDqaJPsi.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-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 iLs" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //
  19 // This class is a part of a package of high level QA monitoring for TRD.
  20 //
  21 // S. Radomski
  22 // radomski@physi.uni-heidelberg.de
  23 // March 2008
  24 //
  25
  26 #include "AliTRDqaJPsi.h"
  27 #include "AliTRDqaAT.h"
  28
  29 #include "TFile.h"
  30 #include "TH1D.h"
  31 #include "TH2D.h"
  32 #include "TChain.h"
  33
  34 #include "AliESDEvent.h"
  35 #include "AliESDtrack.h"
  36 #include "AliKFParticle.h"
  37
  38 #include "TLorentzVector.h"
  39
  40 //const AliTRDqaJPsi::knSteps = 5;
  41
  42 //______________________________________________________________________________
  43
  44 AliTRDqaJPsi::AliTRDqaJPsi()
  45   : AliAnalysisTask("",""),
  46     fChain(0),
  47     fESD(0),
  48     fOutputContainer(0),
  49     fnKFtracks(0)
  50 {
  51   //
  52   // default dummy constructor
  53   //
  54
  55 }
  56 //______________________________________________________________________________
  57
  58 AliTRDqaJPsi:: AliTRDqaJPsi(AliTRDqaJPsi& /*trd*/)
  59   : AliAnalysisTask("",""),
  60     fChain(0),
  61     fESD(0),
  62     fOutputContainer(0),
  63     fnKFtracks(0)
  64 {
  65   //
  66   // Dummy copy constructor
  67   //
  68
  69   //return *this;
  70 }
  71
  72
  73 //______________________________________________________________________________
  74 AliTRDqaJPsi::AliTRDqaJPsi(const char *name)
  75   : AliAnalysisTask(name,""),
  76     fChain(0),
  77     fESD(0),
  78     fOutputContainer(0),
  79     fnKFtracks(0)
  80
  81 {
  82   // Constructor.
  83   // Input slot #0 works with an Ntuple
  84   DefineInput(0, TChain::Class());
  85   // Output slot #0 writes into a TH1 container
  86   DefineOutput(0,  TObjArray::Class()) ;
  87 }
  88
  89 //______________________________________________________________________________
  90 void AliTRDqaJPsi::ConnectInputData(const Option_t *)
  91 {
  92   // Initialisation of branch container and histograms
  93
  94   //AliInfo(Form("*** Initialization of %s", GetName())) ;
  95
  96   fChain = (TChain*)GetInputData(0);
  97   fESD = new AliESDEvent();
  98   fESD->ReadFromTree(fChain);
  99 }
 100
 101 //________________________________________________________________________
 102 void AliTRDqaJPsi::CreateOutputObjects()
 103 {
 104
 105   Int_t c = 0;
 106   fOutputContainer = new TObjArray(100);
 107
 108   const char *charge[2] = {"Neg", "Pos"};
 109
 110   // build histograms
 111
 112   for(Int_t i=0; i<knSteps; i++) {
 113
 114     fStatus[i] = new TH1D(Form("status_%d", i), "status", 32, -0.5, 31.5);
 115     fOutputContainer->AddAt(fStatus[i], c++);
 116
 117     fInvMass[i] = new TH1D(Form("mass_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
 118     fOutputContainer->AddAt(fInvMass[i], c++);
 119
 120     fInvMassVec[i] = new TH1D(Form("massVec_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
 121     fOutputContainer->AddAt(fInvMassVec[i], c++);
 122
 123     fInvMassDiff[i] = new TH1D(Form("massDiff_%d", i), ";m_{inv} (GeV);", 100, -1, 1);
 124     fOutputContainer->AddAt(fInvMassDiff[i], c++);
 125
 126     fAngleSM[i] = new TH1D(Form("angleSM_%d", i), ";#delta SM", 19, -0.5, 18.5);
 127     fOutputContainer->AddAt(fAngleSM[i], c++);
 128
 129     fPtAngle[i] = new TH2D(Form("ptAngle_%d", i), ";p_{T} (GeV/c);#delta SM",
 130                            20, 0, 5, 10, -0.5, 9.5);
 131     fOutputContainer->AddAt(fPtAngle[i], c++);
 132
 133
 134     for(Int_t j=0; j<2; j++) {
 135       fnTracks[j*knSteps+i] =
 136         new TH1D(Form("nTracks%s_%d", charge[j],i), Form("%s;number of tracks",charge[j]), 100, -0.5, 99.5);
 137       fPt[j*knSteps+i] = new TH1D(Form("pt%s_%d", charge[j], i), Form("%s;p_{T} (GeV/c)", charge[j]), 100, 0, 5);
 138       fPID[j*knSteps+i] = new TH1D(Form("pid%s_%d", charge[j], i), ";electron LQ", 100, 0, 1);
 139
 140       fOutputContainer->AddAt(fnTracks[j*knSteps+i], c++);
 141       fOutputContainer->AddAt(fPt[j*knSteps+i], c++);
 142       fOutputContainer->AddAt(fPID[j*knSteps+i], c++);
 143     }
 144   }
 145
 146   //TH2D *fnGoodTracks;
 147
 148   printf("n hist = %d\n", c);
 149 }
 150 //______________________________________________________________________________
 151 void AliTRDqaJPsi::Exec(Option_t *)
 152 {
 153   /*
 154     Selection steps:
 155     - Parameters In and Out
 156     - TRDrefit bit
 157     - TRDpid and quality
 158     - pt > 0.8 GeV/c
 159     - PID > 0.9
 160    */
 161
 162
 163   // Process one event
 164   Long64_t entry = fChain->GetReadEntry() ;
 165   if (!(entry%100)) Info("Exec", "Entry = %ld", entry);
 166
 167   // Processing of one event
 168
 169   if (!fESD) {
 170     //AliError("fESD is not connected to the input!") ;
 171     return ;
 172   }
 173
 174
 175   Int_t nTracks = fESD->GetNumberOfTracks();
 176   Int_t cTracks[2*knSteps] = {0,0,0,0,0,0,0,0,0,0};
 177   fnKFtracks = 0;
 178
 179   // track loop
 180   for(Int_t i=0; i<nTracks; i++) {
 181
 182     //
 183     // track selection
 184     //
 185     // param in and Out
 186     // TRDrefit and TRDPid bit
 187     //
 188
 189     AliESDtrack *track = fESD->GetTrack(i);
 190     const AliExternalTrackParam *paramOut = track->GetOuterParam();
 191     const AliExternalTrackParam *paramIn = track->GetInnerParam();
 192
 193     // long track ..
 194     if (!paramIn) continue;
 195     if (!paramOut) continue;
 196
 197     Int_t step    = 0;
 198     Int_t charge  = (track->Charge() > 0) ? 1 : 0;
 199     UInt_t status = track->GetStatus();
 200     Double_t pt   = track->Pt();
 201     Double_t pid  = track->GetTRDpid(AliPID::kElectron);
 202
 203     Double_t esdPid[5];
 204     track->GetESDpid(esdPid);
 205
 206     // create a kalman particle
 207     Int_t pdg = (charge == 0)? -11 : 11;
 208     for(Int_t k=0; k<knSteps; k++) fInSample[fnKFtracks][k] = 0;
 209
 210     vec[fnKFtracks] = CreateVector(track);
 211     tracks[fnKFtracks] = new AliKFParticle(*track, pdg);
 212     fSM[fnKFtracks] = AliTRDqaAT::GetSector(paramOut->GetAlpha());
 213     fnKFtracks++;
 214
 215     //AliTRDqaAT::PrintPID(track);
 216
 217     // apply the cuts
 218
 219     cTracks[knSteps *charge + step]++;
 220     FillHist(track, step++);
 221
 222     if (!(status & AliESDtrack::kTRDrefit)) continue;
 223
 224     cTracks[knSteps *charge + step]++;
 225     FillHist(track, step++);
 226
 227     if (!(status & AliESDtrack::kTRDpid)) continue;
 228     if (track->GetTRDpidQuality() < 6) continue;
 229
 230     cTracks[knSteps *charge + step]++;
 231     FillHist(track, step++);
 232
 233     if (pt < 0.8) continue;
 234
 235     cTracks[knSteps *charge + step]++;
 236     FillHist(track, step++);
 237
 238     if (pid < 0.3) continue; //
 239     //if (esdPid[AliPID::kElectron] < 0.5) continue;
 240
 241     cTracks[knSteps *charge + step]++;
 242     FillHist(track, step);
 243
 244     for(Int_t k=0; k<2*knSteps; k++) fnTracks[k]->Fill(cTracks[k]);
 245   }
 246
 247   // calculate invariant mass
 248
 249   for(Int_t k=0; k<knSteps; k++) {
 250     for(Int_t i=0; i<fnKFtracks; i++) {
 251       if (!fInSample[i][k]) continue;
 252       for(Int_t j=i+1; j<fnKFtracks; j++) {
 253         if (!fInSample[j][k]) continue;
 254         AliKFParticle jpsi(*(tracks[i]), *(tracks[j]));
 255         TLorentzVector jpsiVec = (*(vec[i])) + (*vec[j]);
 256         fInvMass[k]->Fill(jpsi.GetMass());
 257         fInvMassVec[k]->Fill(jpsiVec.M());
 258         fInvMassDiff[k]->Fill(jpsiVec.M() - jpsi.GetMass());
 259
 260         if (jpsi.GetMass() > 2.5 && jpsi.GetMass() < 3.5) {
 261           Int_t dSM = TMath::Abs(fSM[i] - fSM[j]);
 262           if (dSM > 9) dSM = (18-dSM);
 263           fAngleSM[k]->Fill(dSM);
 264           fPtAngle[k]->Fill(TMath::Hypot(jpsi.GetPx(), jpsi.GetPy()), dSM);
 265         }
 266
 267       }
 268     }
 269   }
 270
 271   PostData(0, fOutputContainer);
 272 }
 273
 274 //______________________________________________________________________________
 275 void AliTRDqaJPsi::Terminate(Option_t *)
 276 {
 277   // save histograms
 278   fOutputContainer = (TObjArray*)GetOutputData(0);
 279
 280   TFile *file = new TFile("outJPsi.root", "RECREATE");
 281   fOutputContainer->Write();
 282
 283   file->Flush();
 284   file->Close();
 285   delete file;
 286
 287   //for(Int_t i=0; i<fOutputContainer->GetEntries(); i++) {
 288   //  TObject *obj = fOu
 289   // }
 290 }
 291
 292 //______________________________________________________________________________
 293
 294 void AliTRDqaJPsi::FillHist(AliESDtrack *track, Int_t step) {
 295
 296   Int_t charge  = (track->Charge() > 0) ? 1 : 0;
 297   UInt_t status = track->GetStatus();
 298   Double_t pt   = track->Pt();
 299   Double_t pid  = track->GetTRDpid(AliPID::kElectron);
 300
 301   Double_t esdPid[5];
 302   track->GetESDpid(esdPid);
 303
 304   Int_t id = charge * knSteps + step;
 305   AliTRDqaAT::FillStatus(fStatus[step], status);
 306   fPt[id]->Fill(pt);
 307   fPID[id]->Fill(pid);
 308   //fPID[id]->Fill(esdPid[AliPID::kElectron]);
 309
 310   fInSample[fnKFtracks-1][step] = 1;
 311 }
 312
 313 //______________________________________________________________________________
 314
 315 TLorentzVector *AliTRDqaJPsi::CreateVector(AliESDtrack *track) {
 316
 317   TLorentzVector *vec = new TLorentzVector();
 318   vec->SetXYZM(track->Px(), track->Py(), track->Pz(), 0);
 319   return vec;
 320 }
 321
 322 //______________________________________________________________________________