PWGCF/EBYE/Fluctuations/AliEbyEHigherMomentsTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Satyajit Jena.                                                 *
   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 //=========================================================================//
  18 //                                                                         //
  19 //           Analysis Task for Net-Charge Higher Moment Analysis           //
  20 //              Author: Satyajit Jena || Nirbhay K. Behera                 //
  21 //                      sjena@cern.ch || nbehera@cern.ch                   //
  22 //                                                                         //
  23 //=========================================================================//
  24
  25 #include "TChain.h"
  26 #include "TList.h"
  27 #include "TFile.h"
  28 #include "TTree.h"
  29 #include "TH1D.h"
  30 #include "TH2F.h"
  31 #include "TCanvas.h"
  32
  33 #include "AliAnalysisTaskSE.h"
  34 #include "AliAnalysisManager.h"
  35
  36 #include "AliESDVertex.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDInputHandler.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODInputHandler.h"
  42 #include "AliESD.h"
  43 #include "AliESDEvent.h"
  44 #include "AliAODEvent.h"
  45 #include "AliStack.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliAODMCHeader.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliMCEvent.h"
  51 #include "AliStack.h"
  52 #include "AliGenHijingEventHeader.h"
  53 #include "AliGenEventHeader.h"
  54 #include "AliPID.h"
  55 #include "AliAODPid.h"
  56 #include "AliPIDResponse.h"
  57 #include "AliAODpidUtil.h"
  58 #include "AliPIDCombined.h"
  59
  60 #include "AliEbyEHigherMomentsTask.h"
  61
  62 using std::cout;
  63 using std::endl;
  64
  65 ClassImp(AliEbyEHigherMomentsTask)
  66
  67
  68 //-----------------------------------------------------------------------
  69 AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
  70 : AliAnalysisTaskSE( name ),
  71   fListOfHistos(0),
  72   fCentralityEstimator("V0M"),
  73   fVxMax(3.),
  74   fVyMax(3.),
  75   fVzMax(15.),
  76   fDCAxy(2.4),
  77   fDCAz(3.),
  78   fPtLowerLimit(0.3),
  79   fPtHigherLimit(1.5),
  80   fEtaLowerLimit(-0.8),
  81   fEtaHigherLimit(0.8),
  82   fTPCNClus(80),
  83   nAODtrackCutBit(128),
  84   fEventCounter(0)
  85 {
  86
  87   for(Int_t i = 0; i < 91; i++)
  88     {
  89       fhNplusNminus[i] = NULL;
  90     }
  91
  92   for ( Int_t i = 0; i < 11; i++) {
  93     fHistQA[i] = NULL;
  94   }
  95
  96   DefineOutput(1, TList::Class()); // Outpput....
  97 }
  98
  99 AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
 100   if(fListOfHistos) delete fListOfHistos;
 101 }
 102
 103 //---------------------------------------------------------------------------------
 104 void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
 105   //For QA-Histograms
 106   fListOfHistos = new TList();
 107   fListOfHistos->SetOwner(kTRUE);
 108   fEventCounter = new TH1D("fEventCounter","EventCounter", 150, 0.5,150.5);
 109   fListOfHistos->Add(fEventCounter);
 110
 111   fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
 112   fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
 113   fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
 114   fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
 115   fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
 116   fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
 117   fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
 118   fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
 119   fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
 120   fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
 121   fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
 122
 123   for(Int_t i = 0; i < 11; i++)
 124     {
 125       fListOfHistos->Add(fHistQA[i]);
 126     }
 127
 128
 129   TString hname;
 130   TString htitle;
 131
 132   for(Int_t i = 0; i < 25; i++)  {
 133     hname  = "hPlusMinusCentBin"; hname+=i;
 134     htitle = "N+ and N- in Cent Bin "; htitle+=i;
 135     fhNplusNminus[i] = new TH2F(hname.Data(),htitle.Data(),1400,0.5,1400.5,1400,0.5,1400.5);
 136     fListOfHistos->Add(fhNplusNminus[i]);
 137   }
 138
 139   for(Int_t i = 25; i < 50; i++) {
 140     hname  = "hPlusMinusCentBin"; hname+=i;
 141     htitle = "N+ and N- in Cent Bin "; htitle+=i;
 142     fhNplusNminus[i] = new TH2F(hname.Data(),htitle.Data(),800,0.5,800.5,800,0.5,800.5);
 143     fListOfHistos->Add(fhNplusNminus[i]);
 144   }
 145
 146   for(Int_t i = 50; i < 91; i++) {
 147     hname  = "hPlusMinusCentBin"; hname+=i;
 148     htitle = "N+ and N- in Cent Bin "; htitle+=i;
 149     fhNplusNminus[i] = new TH2F(hname.Data(),htitle.Data(),500,0.5,500.5,500,0.5,500.5);
 150     fListOfHistos->Add(fhNplusNminus[i]);
 151   }
 152
 153   PostData(1, fListOfHistos);
 154 }
 155
 156
 157 //----------------------------------------------------------------------------------
 158 void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
 159   fEventCounter->Fill(1);
 160   Double_t nPlusCharge = 0.;
 161   Double_t nMinusCharge = 0.;
 162
 163   AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 164   if (!gAOD) {
 165     cout<< "ERROR 01: AOD not found " <<endl;
 166     return;
 167   }
 168
 169
 170   if(fAnalysisType == "AOD") {
 171     AliAODHeader *aodHeader = gAOD->GetHeader();
 172     Double_t cent =  aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 173     Double_t count = (Int_t)cent + 20;
 174     fEventCounter->Fill(count);
 175
 176     if(cent < 0 || cent >= 99.99) return;
 177
 178       fEventCounter->Fill(2);
 179
 180       Int_t nCentrality = (Int_t)cent;
 181
 182       const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
 183
 184       if(vertex) {
 185         Double32_t fCov[6];
 186         vertex->GetCovarianceMatrix(fCov);
 187         if(vertex->GetNContributors() > 0){
 188           if(fCov[5] != 0) {
 189
 190             Double_t lvx = vertex->GetX();
 191             Double_t lvy = vertex->GetY();
 192             Double_t lvz = vertex->GetZ();
 193
 194             fEventCounter->Fill(5);
 195
 196             fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
 197
 198             if(TMath::Abs(lvx) < fVxMax) {
 199               if(TMath::Abs(lvy) < fVyMax) {
 200                 if(TMath::Abs(lvz) < fVzMax) {
 201
 202                   fEventCounter->Fill(6);
 203                   fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
 204
 205                   Int_t ntracks = gAOD->GetNumberOfTracks();
 206
 207                   for(Int_t i = 0; i < ntracks; i++) {
 208                     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(i));
 209                     if(!aodTrack) {
 210                       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 211                       continue;
 212                     }
 213
 214                     if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;
 215
 216                     Float_t dxy, dz;
 217
 218                     dxy = aodTrack->DCA();
 219                     dz  = aodTrack->ZAtDCA();
 220
 221                     Double_t pt = aodTrack->Pt();
 222                     Double_t eta = aodTrack->Eta();
 223                     Double_t nclus = aodTrack->GetTPCNcls();
 224
 225                     if( dxy > fDCAxy ) continue;
 226                     if( dz > fDCAz ) continue;
 227                     if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 228                     if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
 229                     if( nclus < fTPCNClus ) continue;
 230
 231                     fHistQA[8]->Fill(pt);
 232                     fHistQA[9]->Fill(eta);
 233                     fHistQA[10]->Fill(aodTrack->Phi());
 234                     fHistQA[6]->Fill(dxy);
 235                     fHistQA[7]->Fill(dz);
 236
 237                     Short_t gCharge = aodTrack->Charge();
 238
 239
 240                     if(gCharge > 0) nPlusCharge += 1.;
 241                     if(gCharge < 0) nMinusCharge += 1.;
 242
 243                   } //--------- Track Loop
 244                   //  cout << "nCentrality "<<nCentrality<<" nPlusCharge "<< nPlusCharge << " nMinusCharge " << nMinusCharge << endl;
 245                   fhNplusNminus[nCentrality]->Fill(nPlusCharge,nMinusCharge);
 246                   fEventCounter->Fill(7);
 247                 }
 248               }
 249             }
 250           }
 251         }
 252       }
 253     }//AOD--analysis-----
 254
 255   else if(fAnalysisType == "MCAOD") {
 256     TClonesArray *arrayMC = (TClonesArray*) gAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 257     if (!arrayMC) {
 258       AliFatal("Error: MC particles branch not found!\n");
 259     }
 260     AliAODMCHeader *mcHdr=0;
 261     mcHdr=(AliAODMCHeader*)gAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 262     if(!mcHdr) {
 263       printf("MC header branch not found!\n");
 264       return;
 265     }
 266
 267     AliAODHeader *aodHeader = gAOD->GetHeader();
 268     AliCentrality* fcentrality = aodHeader->GetCentralityP();
 269     Double_t cent =fcentrality ->GetCentralityPercentile("V0M");
 270
 271     Int_t nCentrality = (Int_t)cent;
 272     if(cent < 0 || cent >= 91) return;
 273
 274     Bool_t ver = kFALSE;
 275     const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
 276     if(vertex) {
 277       Double32_t fCov[6];
 278       vertex->GetCovarianceMatrix(fCov);
 279       if(vertex->GetNContributors() > 0) {
 280         if(fCov[5] != 0) {
 281
 282             if(TMath::Abs(vertex->GetX()) < fVxMax) {
 283               if(TMath::Abs(vertex->GetY()) < fVyMax) {
 284                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
 285                   ver = kTRUE;
 286                 }
 287               }
 288             }
 289         }
 290         }
 291       }
 292
 293     if(!ver) return;
 294     Int_t nMC = arrayMC->GetEntries();
 295     for (Int_t iMC = 0; iMC < nMC; iMC++) {
 296       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 297       if(!partMC->IsPhysicalPrimary())continue;
 298       if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
 299       if (partMC->Pt() < fPtLowerLimit ||  partMC->Pt() > fPtHigherLimit) continue;
 300
 301       if(partMC->Charge() > 0) nPlusCharge += 1.;
 302       if(partMC->Charge() < 0) nMinusCharge += 1.;
 303
 304     }//MC Track loop-----
 305     fhNplusNminus[nCentrality]->Fill(nPlusCharge,nMinusCharge);
 306     //cout << "Cent=" << nCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl;
 307
 308     }//MCAOD--analysis----
 309   else return;
 310
 311
 312
 313
 314   fEventCounter->Fill(8);
 315   PostData(1, fListOfHistos);
 316
 317 }
 318
 319
 320 //------------------------------------------------------------------------
 321 void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
 322
 323   Info("AliEbyEHigherMomentTask"," Task Successfully finished");
 324
 325 }
 326
 327 //------------------------------------------------------------------------