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