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