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