PWG1/AliAnaFwdDetsQA.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 is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //------------------------------
  17 // Analysis task for quality-assurance
  18 // of forward detectors ESD
  19 //
  20 // 12/06/2009 cvetan.cheshkov@cern.ch
  21 //------------------------------
  22
  23 #include "TChain.h"
  24 #include "TROOT.h"
  25 #include "TFile.h"
  26 #include "TH1F.h"
  27 #include "TH2F.h"
  28 #include "TF1.h"
  29 #include "TCanvas.h"
  30 #include "TVector3.h"
  31 #include "TParticle.h"
  32 #include "AliESDEvent.h"
  33 #include "AliESDv0.h"
  34 #include "AliESDcascade.h"
  35 #include "AliESDMuonTrack.h"
  36 #include "AliESDCaloCluster.h"
  37 #include "AliRun.h"
  38 #include "AliMCEvent.h"
  39 #include "AliStack.h"
  40 #include "AliGenEventHeader.h"
  41 #include "AliPID.h"
  42 #include "AliESDVZERO.h"
  43
  44 #include "AliAnaFwdDetsQA.h"
  45
  46 ClassImp(AliAnaFwdDetsQA)
  47
  48 AliAnaFwdDetsQA::AliAnaFwdDetsQA():
  49 AliAnalysisTaskSE("AliAnaFwdDetsQA"),
  50   fListOfHistos(0),
  51   fT0vtxRec(0),
  52   fT0vtxRecGen(0),
  53   fT0time(0),
  54   fT0time2(0),
  55   fT0mult(0),
  56   fT0vtxRes(0),
  57   fT0ampl(0),
  58   fV0a(0),
  59   fV0c(0),
  60   fV0multA(0),
  61   fV0multC(0),
  62   fV0multAcorr(0),
  63   fV0multCcorr(0),
  64   fV0Acorr(0),
  65   fV0Ccorr(0),
  66   fV0ampl(0)
  67 {
  68   // Default constructor
  69   // Define input and output slots here
  70   // Input slot #0 works with a TChain
  71   DefineInput(0, TChain::Class());
  72   // Output slot #1 TList
  73   DefineOutput(1, TList::Class());
  74 }
  75
  76 AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
  77 AliAnalysisTaskSE(name),
  78   fListOfHistos(0),
  79   fT0vtxRec(0),
  80   fT0vtxRecGen(0),
  81   fT0time(0),
  82   fT0time2(0),
  83   fT0mult(0),
  84   fT0vtxRes(0),
  85   fT0ampl(0),
  86   fV0a(0),
  87   fV0c(0),
  88   fV0multA(0),
  89   fV0multC(0),
  90   fV0multAcorr(0),
  91   fV0multCcorr(0),
  92   fV0Acorr(0),
  93   fV0Ccorr(0),
  94   fV0ampl(0)
  95 {
  96   // Constructor
  97   AliInfo("Constructor AliAnaFwdDetsQA");
  98   // Define input and output slots here
  99   // Input slot #0 works with a TChain
 100   DefineInput(0, TChain::Class());
 101   // Output slot #1 TList
 102   DefineOutput(1, TList::Class());
 103 }
 104
 105 TH1F * AliAnaFwdDetsQA::CreateHisto(const char* name, const char* title,Int_t nBins,
 106                                             Double_t xMin, Double_t xMax,
 107                                             const char* xLabel, const char* yLabel)
 108 {
 109   // create a histogram
 110   TH1F* result = new TH1F(name, title, nBins, xMin, xMax);
 111   result->SetOption("E");
 112   if (xLabel) result->GetXaxis()->SetTitle(xLabel);
 113   if (yLabel) result->GetYaxis()->SetTitle(yLabel);
 114   result->SetMarkerStyle(kFullCircle);
 115   return result;
 116 }
 117
 118 TH1F *AliAnaFwdDetsQA::CreateEffHisto(const TH1F* hGen, const TH1F* hRec)
 119 {
 120   // create an efficiency histogram
 121   Int_t nBins = hGen->GetNbinsX();
 122   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
 123   hEff->SetTitle("");
 124   hEff->SetStats(kFALSE);
 125   hEff->SetMinimum(0.);
 126   hEff->SetMaximum(110.);
 127   hEff->GetYaxis()->SetTitle("#epsilon [%]");
 128
 129   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
 130     Double_t nGenEff = hGen->GetBinContent(iBin);
 131     Double_t nRecEff = hRec->GetBinContent(iBin);
 132     if (nGenEff > 0) {
 133       Double_t eff = nRecEff/nGenEff;
 134       hEff->SetBinContent(iBin, 100. * eff);
 135       Double_t error = sqrt(eff*(1.-eff) / nGenEff);
 136       if (error == 0) error = 0.0001;
 137       hEff->SetBinError(iBin, 100. * error);
 138     }
 139     else {
 140       hEff->SetBinContent(iBin, -100.);
 141       hEff->SetBinError(iBin, 0);
 142     }
 143   }
 144   return hEff;
 145 }
 146
 147
 148 Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
 149 {
 150   // fit a gaussian to a histogram
 151   static TF1* fitFunc = new TF1("fitFunc", "gaus");
 152   fitFunc->SetLineWidth(2);
 153   fitFunc->SetFillStyle(0);
 154   Double_t maxFitRange = 2;
 155
 156   if (histo->Integral() > 50) {
 157     Float_t mean = histo->GetMean();
 158     Float_t rms = histo->GetRMS();
 159     fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms);
 160     fitFunc->SetParameters(mean, rms);
 161     histo->Fit(fitFunc, "QRI0");
 162     histo->GetFunction("fitFunc")->ResetBit(1<<9);
 163     res = TMath::Abs(fitFunc->GetParameter(2));
 164     resError = TMath::Abs(fitFunc->GetParError(2));
 165     return kTRUE;
 166   }
 167   return kFALSE;
 168 }
 169
 170 void AliAnaFwdDetsQA::UserCreateOutputObjects()
 171 {
 172   // Create histograms
 173   AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
 174   // Create output container
 175   fListOfHistos = new TList();
 176
 177   fT0vtxRec = CreateHisto("hT0vtxRec", "Z vertex reconstructed with T0", 100, -25, 25, "Z_{vtx} [cm]", "");
 178   fT0time   = CreateHisto("hT0time", "Time0 reconstruction with T0", 5000, 10000, 20000, "t_{0} [ps]", "");
 179   fT0time2  = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", "");
 180   fT0mult   = CreateHisto("hT0mult", "Total reconstructed multiplicity in T0", 100, -0.5, 99.5, "Multiplicity", "");
 181   fT0vtxRes = CreateHisto("hT0vtxRes", "T0 Z vertex resolution", 100, -25, 25, "Delta(Z_{vtx}) [cm]", "");
 182   fT0ampl   = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5);
 183
 184   fT0vtxRecGen = new TH2F("hT0vtxRecGen", "T0 reconstructed vs generated Z vertex", 100, -25, 25, 100, -25, 25);
 185   fT0vtxRecGen->GetXaxis()->SetTitle("Generated Z vertex");
 186   fT0vtxRecGen->GetYaxis()->SetTitle("Reconstructed Z vertex");
 187   fT0vtxRecGen->SetMarkerStyle(kFullCircle);
 188   fT0vtxRecGen->SetMarkerSize(0.4);
 189
 190   fV0a = CreateHisto("hV0a","Number of fired PMTs (V0A)",65,-0.5,64.5);
 191   fV0c = CreateHisto("hV0c","Number of fired PMTs (V0C)",65,-0.5,64.5);
 192   fV0multA = CreateHisto("hV0multA","Total reconstructed multiplicity (V0A)",100,0.,1000.);
 193   fV0multC = CreateHisto("hV0multC","Total reconstructed multiplicity (V0C)",100,0.,1000.);
 194   fV0ampl  = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5);
 195
 196   fV0multAcorr = new TH2F("hV0multAcorr", "Reconstructed vs generated (primaries only) multiplicity (V0A)",100,0.,500.,100,0.,1000.);
 197   fV0multAcorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
 198   fV0multAcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0A");
 199   fV0multCcorr = new TH2F("hV0multCcorr", "Reconstructed vs generated (primaries only) multiplicity (V0C)",100,0.,500.,100,0.,1000.);
 200   fV0multCcorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
 201   fV0multCcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0C");
 202
 203   fV0Acorr = new TH2F("hV0Acorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0A)",100,0.,500.,65,-0.5,64.5);
 204   fV0Acorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
 205   fV0Acorr->GetYaxis()->SetTitle("Number of fired PMTs in V0A");
 206   fV0Ccorr = new TH2F("hV0Ccorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0C)",100,0.,500.,65,-0.5,64.5);
 207   fV0Ccorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
 208   fV0Ccorr->GetYaxis()->SetTitle("Number of fired PMTs in V0C");
 209
 210   fListOfHistos->Add(fT0vtxRec);
 211   fListOfHistos->Add(fT0time);
 212   fListOfHistos->Add(fT0mult);
 213   fListOfHistos->Add(fT0vtxRecGen);
 214   fListOfHistos->Add(fT0vtxRes);
 215   fListOfHistos->Add(fV0a);
 216   fListOfHistos->Add(fV0c);
 217   fListOfHistos->Add(fV0multA);
 218   fListOfHistos->Add(fV0multC);
 219   fListOfHistos->Add(fV0multAcorr);
 220   fListOfHistos->Add(fV0multCcorr);
 221   fListOfHistos->Add(fV0Acorr);
 222   fListOfHistos->Add(fV0Ccorr);
 223   fListOfHistos->Add(fT0ampl);
 224   fListOfHistos->Add(fV0ampl);
 225   fListOfHistos->Add(fT0time2);
 226 }
 227
 228 void AliAnaFwdDetsQA::UserExec(Option_t */*option*/)
 229 {
 230
 231   AliMCEvent* mcEvent = MCEvent();
 232   if (!mcEvent) {
 233     Printf("ERROR: Could not retrieve MC event");
 234     return;
 235   }
 236
 237   //Primary vertex needed
 238   TArrayF mcvtx(3);
 239   mcEvent->GenEventHeader()->PrimaryVertex(mcvtx);
 240
 241   AliStack *stack = mcEvent->Stack();
 242   if (!stack) {
 243     Printf("ERROR: Could not retrieve MC stack");
 244     return;
 245   }
 246
 247   Int_t nV0A = 0;
 248   Int_t nV0C = 0;
 249   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again
 250     if (!stack->IsPhysicalPrimary(iTracks)) continue;
 251     AliMCParticle* track = mcEvent->GetTrack(iTracks);
 252     TParticle* particle = track->Particle();
 253     if (!particle) continue;
 254     if (track->Charge() == 0) continue;
 255     Double_t eta = particle->Eta();
 256     if (eta > 2.8 && eta < 5.1) {
 257       nV0A++;
 258     }
 259     if (eta > -3.7 && eta < -1.7) {
 260       nV0C++;
 261     }
 262   }
 263
 264   // ESD information
 265   AliVEvent* event = InputEvent();
 266   if (!event) {
 267     Printf("ERROR: Could not retrieve event");
 268     return;
 269   }
 270
 271   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 272   const AliESDTZERO* esdT0 = esd->GetESDTZERO();
 273   Double_t t0zvtx = esdT0->GetT0zVertex();
 274   Double_t t0time = esdT0->GetT0();
 275
 276   fT0vtxRec->Fill(t0zvtx);
 277   fT0time->Fill(t0time);
 278   fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
 279   t0zvtx *= -1.0;
 280   fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
 281
 282   Double_t t0mult = 0;
 283   for(Int_t i = 0; i < 24; i++) {
 284     t0mult += esdT0->GetT0amplitude()[i];
 285     fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
 286   }
 287
 288   fT0mult->Fill(t0mult);
 289   if (t0mult > 10)
 290     fT0time2->Fill(t0time);
 291
 292   AliESDVZERO* esdV0 = esd->GetVZEROData();
 293   fV0a->Fill(esdV0->GetNbPMV0A());
 294   fV0c->Fill(esdV0->GetNbPMV0C());
 295   fV0multA->Fill(esdV0->GetMTotV0A());
 296   fV0multC->Fill(esdV0->GetMTotV0C());
 297
 298   fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
 299   fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
 300   fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
 301   fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
 302
 303   for(Int_t i = 0; i < 64; i++) {
 304     fV0ampl->Fill(esdV0->GetMultiplicity(i));
 305   }
 306   // Post output data.
 307   PostData(1, fListOfHistos);
 308 }
 309
 310 void AliAnaFwdDetsQA::Terminate(Option_t *)
 311 {
 312   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
 313   if (!fListOfHistos) {
 314     Printf("ERROR: fListOfHistos not available");
 315     return;
 316   }
 317
 318   fT0vtxRec = dynamic_cast<TH1F*>(fListOfHistos->At(0));
 319   fT0time = dynamic_cast<TH1F*>(fListOfHistos->At(1));
 320   fT0mult = dynamic_cast<TH1F*>(fListOfHistos->At(2));
 321   fT0vtxRecGen = dynamic_cast<TH2F*>(fListOfHistos->At(3));
 322   fT0vtxRes = dynamic_cast<TH1F*>(fListOfHistos->At(4));
 323
 324   fV0a = dynamic_cast<TH1F*>(fListOfHistos->At(5));
 325   fV0c = dynamic_cast<TH1F*>(fListOfHistos->At(6));
 326   fV0multA = dynamic_cast<TH1F*>(fListOfHistos->At(7));
 327   fV0multC = dynamic_cast<TH1F*>(fListOfHistos->At(8));
 328   fV0multAcorr = dynamic_cast<TH2F*>(fListOfHistos->At(9));
 329   fV0multCcorr = dynamic_cast<TH2F*>(fListOfHistos->At(10));
 330   fV0Acorr = dynamic_cast<TH2F*>(fListOfHistos->At(11));
 331   fV0Ccorr = dynamic_cast<TH2F*>(fListOfHistos->At(12));
 332
 333   fT0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(13));
 334   fV0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(14));
 335   fT0time2 = dynamic_cast<TH1F*>(fListOfHistos->At(15));
 336
 337
 338   // draw the histograms if not in batch mode
 339   if (!gROOT->IsBatch()) {
 340     new TCanvas;
 341     fT0vtxRec->DrawCopy("E");
 342     new TCanvas;
 343     fT0time->DrawCopy("E");
 344     new TCanvas;
 345     fT0mult->DrawCopy("E");
 346     new TCanvas;
 347     fT0vtxRes->DrawCopy("E");
 348     new TCanvas;
 349     fT0vtxRecGen->DrawCopy();
 350     new TCanvas;
 351     fV0a->DrawCopy("E");
 352     new TCanvas;
 353     fV0c->DrawCopy("E");
 354     new TCanvas;
 355     fV0multA->DrawCopy("E");
 356     new TCanvas;
 357     fV0multC->DrawCopy("E");
 358     new TCanvas;
 359     fV0multAcorr->DrawCopy();
 360     new TCanvas;
 361     fV0multCcorr->DrawCopy();
 362     new TCanvas;
 363     fV0Acorr->DrawCopy();
 364     new TCanvas;
 365     fV0Ccorr->DrawCopy();
 366     new TCanvas;
 367     fT0ampl->DrawCopy("E");
 368     new TCanvas;
 369     fV0ampl->DrawCopy("E");
 370     new TCanvas;
 371     fT0time2->DrawCopy("E");
 372   }
 373
 374   // write the output histograms to a file
 375   TFile* outputFile = TFile::Open("FwdDetsQA.root", "recreate");
 376   if (!outputFile || !outputFile->IsOpen())
 377     {
 378       Error("AliAnaFwdDetsQA", "opening output file FwdDetsQA.root failed");
 379       return;
 380     }
 381   fT0vtxRec->Write();
 382   fT0time->Write();
 383   fT0mult->Write();
 384   fT0vtxRecGen->Write();
 385   fT0vtxRes->Write();
 386   fT0ampl->Write();
 387   fT0time2->Write();
 388
 389   fV0a->Write();
 390   fV0c->Write();
 391   fV0multA->Write();
 392   fV0multC->Write();
 393   fV0multAcorr->Write();
 394   fV0multCcorr->Write();
 395   fV0Acorr->Write();
 396   fV0Ccorr->Write();
 397   fV0ampl->Write();
 398
 399   outputFile->Close();
 400   delete outputFile;
 401
 402   //delete esd;
 403   Info("AliAnaFwdDetsQA", "Successfully finished");
 404 }
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466