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