PWGPP/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   if (!esd) {
 279     Printf("ERROR: Could not retrieve esd");
 280     return;
 281   }
 282   const AliESDTZERO* esdT0 = esd->GetESDTZERO();
 283   Double_t t0zvtx = esdT0->GetT0zVertex();
 284   Double_t t0time = esdT0->GetT0();
 285
 286   fT0vtxRec->Fill(t0zvtx);
 287   fT0time->Fill(t0time);
 288   fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
 289   t0zvtx *= -1.0;
 290   fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
 291
 292   Double_t t0mult = 0;
 293   for(Int_t i = 0; i < 24; i++) {
 294     t0mult += esdT0->GetT0amplitude()[i];
 295     fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
 296   }
 297
 298   fT0mult->Fill(t0mult);
 299   if (t0mult > 10)
 300     fT0time2->Fill(t0time);
 301
 302   AliESDVZERO* esdV0 = esd->GetVZEROData();
 303   fV0a->Fill(esdV0->GetNbPMV0A());
 304   fV0c->Fill(esdV0->GetNbPMV0C());
 305   fV0multA->Fill(esdV0->GetMTotV0A());
 306   fV0multC->Fill(esdV0->GetMTotV0C());
 307
 308   fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
 309   fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
 310   fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
 311   fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
 312
 313   for(Int_t i = 0; i < 64; i++) {
 314     fV0ampl->Fill(esdV0->GetMultiplicity(i));
 315   }
 316   // Post output data.
 317   PostData(1, fListOfHistos);
 318 }
 319
 320 void AliAnaFwdDetsQA::Terminate(Option_t *)
 321 {
 322   // Terminate
 323   // Store the output histos
 324   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
 325   if (!fListOfHistos) {
 326     Printf("ERROR: fListOfHistos not available");
 327     return;
 328   }
 329
 330   //delete esd;
 331   Info("AliAnaFwdDetsQA", "Successfully finished");
 332 }
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394