PWGHF/hfe/AliHFEVZEROEventPlane.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 // VZERO event plane task for 2010
  18 // Gain equalization + Recentering
  19 // Need a root file
  20 //
  21 // Author:
  22 //   Raphaelle Bailhache <R.Bailhache@gsi.de>
  23 //
  24 //
  25
  26 #include "AliHFEVZEROEventPlane.h"
  27
  28 // AliRoot includes
  29 #include "AliAnalysisManager.h"
  30 #include "AliInputEventHandler.h"
  31 #include "AliVEvent.h"
  32 //#include "AliESDEvent.h"
  33 #include "AliCentrality.h"
  34 //#include "AliESDVZERO.h"
  35 #include "AliVVZERO.h"
  36 #include "TFile.h"
  37 #include "AliOADBContainer.h"
  38 #include "TH2F.h"
  39 #include "TF1.h"
  40 #include "TList.h"
  41 #include "TChain.h"
  42 #include "THnSparse.h"
  43 #include "TString.h"
  44 #include "TVector2.h"
  45 #include "AliEventplane.h"
  46
  47
  48 // STL includes
  49 #include <iostream>
  50 using namespace std;
  51
  52
  53 //_____________________________________________________________________________
  54 AliHFEVZEROEventPlane::AliHFEVZEROEventPlane():
  55   TNamed(),
  56   fEventPlaneV0A(-100.0),
  57   fEventPlaneV0C(-100.0),
  58   fEventPlaneV0(-100.0),
  59   fESD(0),
  60   fRun(-1),
  61   fMultV0(0x0),
  62   fV0Cpol(100),
  63   fV0Apol(100),
  64   fnamefile(TString("")),
  65   fOutputList(0x0),
  66   fMultV0Before(0x0),
  67   fMultV0After(0x0)
  68 {
  69   //
  70   // Default constructor (should not be used)
  71   //
  72
  73   for(Int_t k = 0; k < fgknCentrBin; k++) {
  74     for(Int_t iside = 0; iside < 2; iside++) {
  75       for(Int_t icoord = 0; icoord < 2; icoord++) {
  76         fQBefore[k][iside][icoord] = 0x0;
  77         fQAfter[k][iside][icoord] = 0x0;
  78       }
  79     }
  80   }
  81
  82 }
  83
  84 //______________________________________________________________________________
  85 AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const char *name, const Char_t *title):
  86   TNamed(name,title),
  87   fEventPlaneV0A(-100.0),
  88   fEventPlaneV0C(-100.0),
  89   fEventPlaneV0(-100.0),
  90   fESD(0),
  91   fRun(-1),
  92   fMultV0(0x0),
  93   fV0Cpol(100),
  94   fV0Apol(100),
  95   fnamefile(TString("")),
  96   fOutputList(0x0),
  97   fMultV0Before(0x0),
  98   fMultV0After(0x0)
  99 {
 100   //
 101   // Constructor
 102   //
 103
 104   TString namelist("QA_hist_");
 105   namelist += name;
 106   fOutputList = new TList();
 107   fOutputList->SetName((const char*)namelist);
 108   fOutputList->SetOwner(kTRUE);
 109
 110   // Multiplicity
 111   fMultV0Before = new TProfile("MultV0Before","",64,0,64);
 112   fMultV0Before->Sumw2();
 113   fMultV0After = new TProfile("MultV0After","",64,0,64);
 114   fMultV0After->Sumw2();
 115   fOutputList->Add(fMultV0Before);
 116   fOutputList->Add(fMultV0After);
 117
 118   // Recentering
 119   for(Int_t k = 0; k < fgknCentrBin; k++) {
 120     for(Int_t iside = 0; iside < 2; iside++) {
 121       for(Int_t icoord = 0; icoord < 2; icoord++) {
 122
 123         TString namecontbefore("");
 124         TString namecontafter("");
 125
 126         if(iside==0 && icoord==0) {
 127           namecontbefore += "hQxc2_";
 128           namecontbefore += k;
 129         }
 130         else if(iside==1 && icoord==0){
 131           namecontbefore += "hQxa2_";
 132           namecontbefore += k;
 133         }
 134         else if(iside==0 && icoord==1){
 135           namecontbefore += "hQyc2_";
 136           namecontbefore += k;
 137         }
 138         else if(iside==1 && icoord==1) {
 139           namecontbefore += "hQya2_";
 140           namecontbefore += k;
 141         }
 142         //
 143         if(iside==0 && icoord==0){
 144           namecontafter += "hQxc2_";
 145           namecontafter += k;
 146           namecontafter += "_after";
 147         }
 148         else if(iside==1 && icoord==0){
 149           namecontafter += "hQxa2_";
 150           namecontafter += k;
 151           namecontafter += "_after";
 152         }
 153         else if(iside==0 && icoord==1) {
 154           namecontafter += "hQyc2_";
 155           namecontafter += k;
 156           namecontafter += "_after";
 157         }
 158         else if(iside==1 && icoord==1) {
 159           namecontafter += "hQya2_";
 160           namecontafter += k;
 161           namecontafter += "_after";
 162         }
 163         //
 164         fQBefore[k][iside][icoord] = new TH1F(((const char*)namecontbefore),"",1600,-800.0,800.0);
 165         fQBefore[k][iside][icoord]->Sumw2();
 166         fOutputList->Add(fQBefore[k][iside][icoord]);
 167         fQAfter[k][iside][icoord] = new TH1F(((const char*)namecontafter),"",1600,-800.0,800.0);
 168         fQAfter[k][iside][icoord]->Sumw2();
 169         fOutputList->Add(fQAfter[k][iside][icoord]);
 170         //
 171       }
 172     }
 173   }
 174
 175 }
 176 //____________________________________________________________
 177 AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const AliHFEVZEROEventPlane &ref):
 178   TNamed(ref),
 179   fEventPlaneV0A(ref.fEventPlaneV0A),
 180   fEventPlaneV0C(ref.fEventPlaneV0C),
 181   fEventPlaneV0(ref.fEventPlaneV0),
 182   fESD(0),
 183   fRun(-1),
 184   fMultV0(0x0),
 185   fV0Cpol(100),
 186   fV0Apol(100),
 187   fnamefile(ref.fnamefile),
 188   fOutputList(0x0),
 189   fMultV0Before(0x0),
 190   fMultV0After(0x0)
 191 {
 192   //
 193   // Copy Constructor
 194   //
 195   ref.Copy(*this);
 196 }
 197
 198 //____________________________________________________________
 199 AliHFEVZEROEventPlane &AliHFEVZEROEventPlane::operator=(const AliHFEVZEROEventPlane &ref){
 200   //
 201   // Assignment operator
 202   //
 203   if(this == &ref)
 204     ref.Copy(*this);
 205   return *this;
 206 }
 207
 208 //____________________________________________________________
 209 void AliHFEVZEROEventPlane::Copy(TObject &o) const {
 210   //
 211   // Copy into object o
 212   //
 213   AliHFEVZEROEventPlane &target = dynamic_cast<AliHFEVZEROEventPlane &>(o);
 214
 215   target.fEventPlaneV0A = fEventPlaneV0A;
 216   target.fEventPlaneV0C = fEventPlaneV0C;
 217   target.fEventPlaneV0 = fEventPlaneV0;
 218   target.fnamefile = fnamefile;
 219
 220 }
 221
 222 //__________________________________________________________________
 223 AliHFEVZEROEventPlane::~AliHFEVZEROEventPlane(){
 224   //
 225   // Destruktor
 226   //
 227   if(fOutputList){
 228     fOutputList->Delete();
 229     delete fOutputList;
 230   }
 231   fOutputList = 0x0;
 232
 233 }
 234 //______________________________________________________________________________
 235 void AliHFEVZEROEventPlane::ProcessEvent(AliVEvent *event)
 236 {
 237   //
 238   // Process the event
 239   //
 240
 241   // Reset
 242   fEventPlaneV0A = -100.0;
 243   fEventPlaneV0C = -100.0;
 244   fEventPlaneV0  = -100.0;
 245
 246   //
 247   Int_t run = event->GetRunNumber();
 248   Bool_t ok = kTRUE;
 249   if(run != fRun){
 250     if(fMultV0) delete fMultV0;
 251     fMultV0 = 0x0;
 252     // Load the calibrations run dependent
 253     if(!OpenInfoCalbration(run)) ok = kFALSE;
 254     fRun=run;
 255   }
 256
 257   if(ok) Analyze(event);
 258
 259 }
 260
 261 //________________________________________________________________________
 262 void AliHFEVZEROEventPlane::Analyze(AliVEvent* esdEvent)
 263 {
 264   //
 265   // Do VZERO calibration + centering
 266   //
 267
 268   if(!fMultV0) {
 269     //printf("Did not find calibration VZERO\n");
 270     return;
 271   }
 272   //printf("HERE!!!\n");
 273
 274   //Centrality
 275   Float_t v0Centr  = -10.;
 276   AliCentrality* centrality = esdEvent->GetCentrality();
 277   if (centrality){
 278     v0Centr  = centrality->GetCentralityPercentile("V0M");
 279   }
 280
 281   // Analyse only for 0-80% PbPb collisions
 282   Int_t iC = -1;
 283   if (v0Centr >0 && v0Centr < 80){
 284     if(v0Centr < 5) iC = 0;
 285     else if(v0Centr < 10) iC = 1;
 286     else if(v0Centr < 20) iC = 2;
 287     else if(v0Centr < 30) iC = 3;
 288     else if(v0Centr < 40) iC = 4;
 289     else if(v0Centr < 50) iC = 5;
 290     else if(v0Centr < 60) iC = 6;
 291     else if(v0Centr < 70) iC = 7;
 292     else iC = 8;
 293
 294     //general
 295     Double_t qxa2 = 0, qya2 = 0;
 296     Double_t qxc2 = 0, qyc2 = 0;
 297
 298     //V0 info
 299     AliVVZERO* esdV0 = esdEvent->GetVZEROData();
 300
 301     for (Int_t iv0 = 0; iv0 < 64; iv0++) {
 302       Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
 303       Float_t multv0 = esdV0->GetMultiplicity(iv0);
 304       if (iv0 < 32){
 305         //printf("Value %f\n",fMultV0->GetBinContent(iv0+1));
 306         if(fMultV0->GetBinContent(iv0+1)>0.0) {
 307           qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
 308           qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
 309           fMultV0After->Fill(iv0,multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1));
 310         }
 311       } else {
 312         if(fMultV0->GetBinContent(iv0+1)>0.0) {
 313           qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
 314           qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
 315           fMultV0After->Fill(iv0,multv0*fV0Apol/fMultV0->GetBinContent(iv0+1));
 316         }
 317       }
 318       fMultV0Before->Fill(iv0,multv0);
 319     }
 320
 321     // Fill histos
 322     fQBefore[iC][1][0] -> Fill(qxa2);
 323     fQBefore[iC][1][1] -> Fill(qya2);
 324     fQBefore[iC][0][0] -> Fill(qxc2);
 325     fQBefore[iC][0][1] -> Fill(qyc2);
 326
 327     //grab for each centrality the proper histo with the Qx and Qy to do the recentering
 328     Double_t qxamean2 = fMeanQ[iC][1][0];
 329     Double_t qxarms2  = fWidthQ[iC][1][0];
 330     Double_t qyamean2 = fMeanQ[iC][1][1];
 331     Double_t qyarms2  = fWidthQ[iC][1][1];
 332
 333     Double_t qxcmean2 = fMeanQ[iC][0][0];
 334     Double_t qxcrms2  = fWidthQ[iC][0][0];
 335     Double_t qycmean2 = fMeanQ[iC][0][1];
 336     Double_t qycrms2  = fWidthQ[iC][0][1];
 337
 338     if((TMath::Abs(qxarms2) < 0.00000001) || (TMath::Abs(qyarms2) < 0.00000001) || (TMath::Abs(qxcrms2) < 0.00000001) || (TMath::Abs(qycrms2) < 0.00000001)) return;
 339
 340     Double_t qxaCor2 = (qxa2 - qxamean2)/qxarms2;
 341     Double_t qyaCor2 = (qya2 - qyamean2)/qyarms2;
 342     Double_t qxcCor2 = (qxc2 - qxcmean2)/qxcrms2;
 343     Double_t qycCor2 = (qyc2 - qycmean2)/qycrms2;
 344
 345     Double_t qxCor2 = qxaCor2 + qxcCor2;
 346     Double_t qyCor2 = qyaCor2 + qycCor2;
 347
 348     // Fill histos
 349     fQAfter[iC][1][0] -> Fill(qxaCor2);
 350     fQAfter[iC][1][1] -> Fill(qyaCor2);
 351     fQAfter[iC][0][0] -> Fill(qxcCor2);
 352     fQAfter[iC][0][1] -> Fill(qycCor2);
 353
 354     Double_t evPlAngV0ACor2 = TVector2::Phi_0_2pi(TMath::ATan2(qyaCor2, qxaCor2))/2.;
 355     Double_t evPlAngV0CCor2 = TVector2::Phi_0_2pi(TMath::ATan2(qycCor2, qxcCor2))/2.;
 356     Double_t evPlAngV0Cor2  = TVector2::Phi_0_2pi(TMath::ATan2(qyCor2, qxCor2))/2.;
 357
 358     fEventPlaneV0A = evPlAngV0ACor2;
 359     fEventPlaneV0C = evPlAngV0CCor2;
 360     fEventPlaneV0  = evPlAngV0Cor2;
 361
 362     //printf("Eventplane V0A %f, V0C %f\n",fEventPlaneV0A,fEventPlaneV0C);
 363
 364   }
 365 }
 366 //_____________________________________________________________________________
 367 Bool_t AliHFEVZEROEventPlane::OpenInfoCalbration(Int_t run)
 368 {
 369   //
 370   // Take the calibration coefficients
 371   //
 372
 373   //printf("Name of the file %s\n",(const char*)fnamefile);
 374
 375   TFile *foadb = TFile::Open(fnamefile.Data());
 376   if(!foadb){
 377     printf("OADB file %s cannot be opened\n",fnamefile.Data());
 378     return kFALSE;
 379   }
 380
 381   //printf("test\n");
 382
 383   AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
 384   if(!cont){
 385     printf("OADB object hMultV0BefCorr is not available in the file\n");
 386     return kFALSE;
 387   }
 388
 389   if(!(cont->GetObject(run))){
 390     printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
 391     return kFALSE;
 392   }
 393   //TProfile *multV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
 394   TProfile *multV0 = ((TProfile *) cont->GetObject(run));
 395   if(!multV0) return kFALSE;
 396   fMultV0 = (TProfile *) multV0->Clone();
 397   fMultV0->SetDirectory(0);
 398
 399   TF1 *fpol0 = new TF1("fpol0","pol0");
 400   fMultV0->Fit(fpol0,"Q0","",0,31);
 401   fV0Cpol = fpol0->GetParameter(0);
 402   fMultV0->Fit(fpol0,"Q0","",32,64);
 403   fV0Apol = fpol0->GetParameter(0);
 404
 405   for(Int_t iside=0;iside<2;iside++){
 406     for(Int_t icoord=0;icoord<2;icoord++){
 407       for(Int_t i=0;i  < fgknCentrBin;i++){
 408
 409         TString namecont("");
 410         if(iside==0 && icoord==0) {
 411           namecont += "hQxc2_";
 412           namecont += i;
 413         }
 414         else if(iside==1 && icoord==0){
 415           namecont += "hQxa2_";
 416           namecont += i;
 417         }
 418         else if(iside==0 && icoord==1){
 419           namecont += "hQyc2_";
 420           namecont += i;
 421         }
 422         else if(iside==1 && icoord==1) {
 423           namecont += "hQya2_";
 424           namecont += i;
 425         }
 426
 427         cont = (AliOADBContainer*) foadb->Get(namecont);
 428         if(!cont){
 429           printf("OADB object %s is not available in the file\n",(const char*)namecont);
 430           delete fpol0;
 431           return kFALSE;
 432         }
 433
 434         if(!(cont->GetObject(run))){
 435           printf("OADB object %s is not available for run %i\n",(const char*)namecont,run);
 436           delete fpol0;
 437           return kFALSE;
 438         }
 439         fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
 440         fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
 441       }
 442     }
 443   }
 444
 445   delete fpol0;
 446   foadb->Close();
 447
 448   return kTRUE;
 449
 450 }