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