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