PWGHF/hfe/AliHFEVZEROEventPlane.cxx

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