PWGCF/EBYE/Fluctuations/AliEbyEFluctuationAnalysisTask.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TCanvas.h"
   6 #include "TParticle.h"
   7 #include "TObjArray.h"
   8
   9 #include "AliAnalysisTask.h"
  10 #include "AliAnalysisManager.h"
  11
  12 #include "AliStack.h"
  13 #include "AliMCEvent.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDInputHandler.h"
  16 #include "AliESDtrackCuts.h"
  17 #include "AliCentrality.h"
  18
  19 #include "AliEbyEFluctuationAnalysisTask.h"
  20
  21 // Event by event charge fluctuation analysis
  22 // Authors: Satyajit Jena and Panos Cristakoglou
  23 //
  24
  25 ClassImp(AliEbyEFluctuationAnalysisTask)
  26
  27 //________________________________________________________________________
  28 AliEbyEFluctuationAnalysisTask::AliEbyEFluctuationAnalysisTask(const char *name)
  29   : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
  30     fHistEventStats(0), fHistCentrality(0),
  31     fHistNMultMC(0), fHistNPlusNMinusMC(0),
  32     fESDtrackCuts(0),
  33     fAnalysisType("ESD"), fAnalysisMode("TPC"),
  34     fCentralityEstimator("V0M"), fCentralityBins20(kFALSE),
  35     fVxMax(3.), fVyMax(3.), fVzMax(10.) {
  36   // Constructor
  37
  38   for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
  39     fHistNMult[iBin-1] = NULL;
  40     fHistNPlusNMinus[iBin-1] = NULL;
  41   }
  42
  43   // Define input and output slots here
  44   // Input slot #0 works with a TChain
  45   DefineInput(0, TChain::Class());
  46   // Output slot #0 id reserved by the base class for AOD
  47   // Output slot #1 writes into a TH1 container
  48   DefineOutput(1, TList::Class());
  49 }
  50
  51 //________________________________________________________________________
  52 void AliEbyEFluctuationAnalysisTask::UserCreateOutputObjects() {
  53   // Create histograms
  54   // Called once
  55
  56   fOutputList = new TList();
  57   fOutputList->SetOwner();
  58
  59   //Event stats.
  60   TString gCutName[4] = {"Total","Offline trigger",
  61                          "Vertex","Analyzed"};
  62   fHistEventStats = new TH1F("fHistEventStats",
  63                              "Event statistics;;N_{events}",
  64                              4,0.5,4.5);
  65   for(Int_t i = 1; i <= 4; i++)
  66     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
  67   fOutputList->Add(fHistEventStats);
  68
  69   //ESD analysis
  70   if(fAnalysisType == "ESD") {
  71     fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
  72                                nCentralityBins,0.5,nCentralityBins+0.5);
  73     fOutputList->Add(fHistCentrality);
  74
  75     TString histName;
  76     for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
  77       histName = "fHistNMult"; histName += "Centrality"; histName += iBin;
  78       fHistNMult[iBin-1] = new TH1F(histName.Data(),
  79                                     ";N_{mult.}",
  80                                     500,0,3000);
  81       fOutputList->Add(fHistNMult[iBin-1]);
  82     }
  83     for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
  84       histName = "fHistNPlusNMinus"; histName += "Centrality"; histName += iBin;
  85       fHistNPlusNMinus[iBin-1] = new TH2F(histName.Data(),
  86                                           ";N_{+};N_{-}",
  87                                           2000,0.5,2000.5,2000,0.5,2000.5);
  88       fOutputList->Add(fHistNPlusNMinus[iBin-1]);
  89     }
  90   }//ESD analysis
  91   else if(fAnalysisType == "MC") {
  92     TString histName = "fHistNMultMC";
  93     fHistNMultMC = new TH1F(histName.Data(),
  94                             ";N_{mult.}",
  95                             600,0,6000);
  96     fOutputList->Add(fHistNMultMC);
  97
  98     histName = "fHistNPlusNMinusMC";
  99     fHistNPlusNMinusMC = new TH2F(histName.Data(),
 100                                   ";N_{+};N_{-}",
 101                                   3000,0.5,3000.5,3000,0.5,3000.5);
 102     fOutputList->Add(fHistNPlusNMinusMC);
 103   }//MC analysis
 104 }
 105
 106 //________________________________________________________________________
 107 void AliEbyEFluctuationAnalysisTask::UserExec(Option_t *) {
 108   // Main loop
 109   // Called for each event
 110
 111   Int_t nPlus = 0, nMinus = 0;
 112
 113   // Post output data.
 114   //ESD analysis
 115   if(fAnalysisType == "ESD") {
 116     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 117     if (!fESD) {
 118       printf("ERROR: fESD not available\n");
 119       return;
 120     }
 121
 122     fHistEventStats->Fill(1); //all events
 123
 124     //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 125     //if(isSelected) {
 126       fHistEventStats->Fill(2); //triggered + centrality
 127       Printf("Event accepted");
 128
 129       //Centrality stuff
 130       AliCentrality *centrality = fESD->GetCentrality();
 131       Int_t nCentrality = 0;
 132       if(fCentralityBins20)
 133         nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());
 134       else
 135         nCentrality = centrality->GetCentralityClass10(fCentralityEstimator.Data());
 136
 137       if((nCentrality < 0)||(nCentrality > 19)) return;
 138
 139       if(fAnalysisMode == "TPC") {
 140         const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
 141         if(vertex) {
 142           if(vertex->GetNContributors() > 0) {
 143             if(vertex->GetZRes() != 0) {
 144               fHistEventStats->Fill(3); //events with a proper vertex
 145               if(TMath::Abs(vertex->GetX()) < fVxMax) {
 146                 if(TMath::Abs(vertex->GetY()) < fVyMax) {
 147                   if(TMath::Abs(vertex->GetZ()) < fVzMax) {
 148                     fHistEventStats->Fill(4); //analyzed events
 149
 150                     //Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
 151                     //centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
 152                     //nCentrality,fESD->GetNumberOfTracks());
 153
 154                     Int_t nAcceptedTracks = 0;
 155                     TObjArray *gTrackArray = 0;
 156                     if(fESDtrackCuts)
 157                       gTrackArray = fESDtrackCuts->GetAcceptedTracks(fESD,kTRUE);
 158                     if(gTrackArray) {
 159                       nAcceptedTracks = gTrackArray->GetEntries();
 160
 161                       AliESDtrack* track = 0;
 162                       // Track loop to fill a pT spectrum
 163                       for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) {
 164                         track = dynamic_cast<AliESDtrack *>(gTrackArray->At(iTracks));
 165                         if (!track) {
 166                           printf("ERROR: Could not receive track %d\n", iTracks);
 167                           continue;
 168                         }
 169
 170                         Short_t gCharge = track->Charge();
 171                         if(gCharge > 0) nPlus += 1;
 172                         if(gCharge < 0) nMinus += 1;
 173                       }//track loop
 174                     }//TObjArray valid object
 175                     //if((nCentrality >= 1)&&(nCentrality <= 20)) {
 176
 177                     fHistCentrality->Fill(nCentrality);
 178                     fHistNPlusNMinus[nCentrality-1]->Fill(nPlus,nMinus);
 179                     fHistNMult[nCentrality-1]->Fill(nPlus+nMinus);
 180                     //}
 181                   }//Vz cut
 182                 }//Vy cut
 183               }//Vx cut
 184             }//Vz resolution
 185           }//number of contributors
 186         }//valid vertex
 187       }//TPC analysis mode
 188       //}//physics selection
 189   }//ESD analysis level
 190
 191   //MC analysis
 192   if(fAnalysisType == "MC") {
 193     AliMCEvent* mcEvent = MCEvent();
 194     if (!mcEvent) {
 195       Printf("ERROR: Could not retrieve MC event");
 196       return;
 197     }
 198     AliStack* stack = mcEvent->Stack();
 199     if (!stack) {
 200       Printf("ERROR: Could not retrieve MC stack");
 201       return;
 202     }
 203
 204     fHistEventStats->Fill(1);
 205     fHistEventStats->Fill(2);
 206     fHistEventStats->Fill(3);
 207     fHistEventStats->Fill(4); //analyzed events
 208     for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
 209       AliVParticle* track = mcEvent->GetTrack(iTracks);
 210       if (!track) {
 211         Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
 212         continue;
 213       }
 214
 215       if(!stack->IsPhysicalPrimary(iTracks)) continue;
 216       if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;
 217       if(TMath::Abs(track->Eta()) > 0.8) continue;
 218
 219       Short_t gCharge = track->Charge();
 220       if(gCharge > 0) nPlus += 1;
 221       if(gCharge < 0) nMinus += 1;
 222     }//particle loop
 223     fHistNPlusNMinusMC->Fill(nPlus,nMinus);
 224     fHistNMultMC->Fill(nPlus+nMinus);
 225
 226   }//MC analysis level
 227
 228
 229   PostData(1, fOutputList);
 230 }
 231
 232 //________________________________________________________________________
 233 void AliEbyEFluctuationAnalysisTask::Terminate(Option_t *) {
 234   // Draw result to the screen
 235   // Called once at the end of the query
 236
 237   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 238   if (!fOutputList) {
 239     printf("ERROR: Output list not available\n");
 240     return;
 241   }
 242 }