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