PWGCF/EBYE/Fluctuations/AliEbyEMultFluctuationTask.cxx

   1
   2 #include "TChain.h"
   3 #include "TTree.h"
   4 #include "TNtuple.h"
   5 #include "TH1F.h"
   6 #include "TH2D.h"
   7 #include "TH1D.h"
   8 #include "TCanvas.h"
   9 #include "TArrayD.h"
  10 #include "AliAnalysisTask.h"
  11 #include "AliAnalysisManager.h"
  12 #include "AliAODHeader.h"
  13 #include "AliAODTrack.h"
  14 #include "AliAODTracklets.h"
  15 #include "AliAODpidUtil.h"
  16 #include "AliVTrack.h"
  17 #include "AliAODEvent.h"
  18 #include "AliAODInputHandler.h"
  19 #include "AliCentrality.h"
  20 #include "AliEbyEMultFluctuationTask.h"
  21 #include "AliAODVertex.h"
  22 #include "AliAODVZERO.h"
  23 #include "AliCentrality.h"
  24 #include "AliPID.h"
  25 #include "AliPIDResponse.h"
  26 #include "AliPIDCombined.h"
  27 #include "AliAODpidUtil.h"
  28 #include "TRandom.h"
  29 #include "AliTPCPIDResponse.h"
  30 #include "AliTRDPIDResponse.h"
  31
  32
  33 ClassImp(AliEbyEMultFluctuationTask)
  34
  35 //________________________________________________________________________
  36 AliEbyEMultFluctuationTask::AliEbyEMultFluctuationTask(const char *name)
  37 : AliAnalysisTaskSE(name), fAOD(0), fOutputList(0),fHistNchPt(0),fHistNchEta(0),fHistNchEtaCent(0),fHistNchPhi(0),fHistDCAxy(0),fHistDCAz(0),fHistnclus(0),fHistchi2ndf(0),fDeDx(NULL),fHistMult(0),fAODVertex(0),fPID(0),fEventCounter(0),histcounter(0),fHistVz(0),fHistchi2ndfvscs(0),My_ntuple(0),fHistMultV0A(0),fHistMultV0C(0),fHistMultV0total(0),fCentralityCounter(0),fCentralityEstimator("V0M"),Counter(0)
  38
  39
  40
  41 {
  42   // Constructor
  43   for(Int_t ibin=0;ibin<91;ibin++)
  44                 {
  45                         fMult[ibin]=NULL;
  46                 }
  47                 for(Int_t jbin=0;jbin<46;jbin++)
  48                 {
  49                 fMultTwo[jbin]=NULL;
  50                 }
  51                 for(Int_t kbin=0;kbin<15;kbin++)
  52                 {
  53                 fMultFive[kbin]=NULL;
  54                 }
  55
  56
  57   DefineInput(0, TChain::Class());
  58
  59   DefineOutput(1, TList::Class());
  60 }
  61
  62 //________________________________________________________________________
  63 Bool_t AliEbyEMultFluctuationTask :: SelectEvent(AliAODVertex* vertex)
  64 {
  65 if(vertex)
  66
  67       if(vertex->GetNContributors() < 0) return kFALSE;
  68
  69
  70           Double_t lvx = vertex->GetX();
  71           Double_t lvy = vertex->GetY();
  72           Double_t lvz = vertex->GetZ();
  73
  74           fEventCounter->Fill(3);
  75
  76
  77           if(TMath::Abs(lvx) > 0.3)  return kFALSE;
  78             if(TMath::Abs(lvy) > 0.3) return kFALSE;
  79             if(TMath::Abs(lvz) > 10) return kFALSE;
  80             if(vertex->GetType()==AliAODVertex::kPileupSPD) return kFALSE;
  81
  82                 fEventCounter->Fill(5);
  83                 fHistVz->Fill(lvz);
  84 return kTRUE;
  85 }
  86 //_____________________________________________________________________________
  87 Int_t AliEbyEMultFluctuationTask :: SelectTrack(AliAODTrack* track)
  88 {
  89
  90         Double_t eta = track->Eta();
  91         Double_t phi = track->Phi();
  92
  93      if(TMath :: Abs(eta)>0.8) return 0;
  94
  95  if(!track->TestFilterBit(128)) return 0;
  96
  97  Float_t dxy, dz;
  98
  99   dxy = track->DCA();
 100  dz  = track->ZAtDCA();
 101                     if(TMath ::Abs(dxy) >2.4 || TMath ::Abs(dz)>3.0) return 0;
 102                   // cout<<dxy<<dz<<endl;
 103
 104
 105   Double_t nclus = track->GetTPCClusterInfo(2,1);
 106   if(nclus<80) return 0;
 107
 108 //chi2cut
 109
 110  Double_t chi2ndf = track->Chi2perNDF();
 111   if(chi2ndf>4) return 0;
 112
 113 return 1;
 114
 115
 116 }
 117
 118 //___________________________________________________________________________________
 119 void AliEbyEMultFluctuationTask::UserCreateOutputObjects()
 120 {
 121   fOutputList = new TList();
 122
 123
 124
 125   My_ntuple = new TNtuple("My_ntuple", "My_ntuple", "Mult:Mult1:Mult2:nCentrality:fV0total:nBin2:nBin5:spdmult0:spdmult1:tracklets:run");
 126     fOutputList->Add(My_ntuple);
 127
 128 fOutputList->SetOwner(kTRUE);
 129 fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
 130 fOutputList->Add(fEventCounter);
 131
 132
 133 histcounter = new TH1D("histcounter","histcounter", 10, 0.5,10.5);
 134 fOutputList->Add(histcounter);
 135
 136
 137 TString BinLebel;
 138 fCentralityCounter = new TH1D("fEventCounter","EventCounter", 100, 0.5,100.5);
 139 TAxis *xaxis = fCentralityCounter->GetXaxis();
 140  for(Int_t ibin=0;ibin<100;ibin++)
 141 {
 142 BinLebel="cent";BinLebel+=ibin;
 143 xaxis->SetBinLabel(ibin,BinLebel.Data());
 144
 145 }
 146 fOutputList->Add(fCentralityCounter);
 147
 148 fHistMultV0A = new TH1F("fHistMultV0A","Mult-VOA",22000,0,22000);
 149 fHistMultV0C = new TH1F("fHistMultV0C","Mult-VOC",22000,0,22000);
 150 fHistMultV0total = new TH1F("fHistMultV0total","Mult-VOTotal",22000,0,22000);
 151
 152
 153  fHistNchEta = new TH1D("fHistNchEta", "#eta distribution of charged tracks(for hybrid)", 200, -0.8, 0.8);
 154 fHistNchEta->GetXaxis()->SetTitle("#eta");
 155 fHistNchEta->GetYaxis()->SetTitle("dNch/d#eta");
 156
 157
 158 fHistNchEtaCent = new TH1D("fHistNchEtaCent", "Central #eta distribution of charged tracks(for hybrid)", 200, -0.8, 0.8);
 159 fHistNchEtaCent->GetXaxis()->SetTitle("#eta");
 160 fHistNchEtaCent->GetYaxis()->SetTitle("dNch/d#eta(Central)");
 161
 162
 163
 164 fHistNchPhi = new TH1D("fHistNchPhi", "#phi distribution of charged tracks(for hybrid)", 200, -10, 10);
 165 fHistNchPhi->GetXaxis()->SetTitle("#phi");
 166 fHistNchPhi->GetYaxis()->SetTitle("dNch/d#phi");
 167
 168 fHistNchPt = new TH1D("fHistNchPt", "P_{T} distribution(charged tracks for hybrid)", 150, 0.1, 3.1);
 169 fHistNchPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 170 fHistNchPt->GetYaxis()->SetTitle("dNch/dP_{T}(c/GeV)");
 171
 172
 173  fHistchi2ndfvscs = new TH2D("fHistchi2ndfvscs","Chi2/ndf from Gaussian fits of Multiplicity Distribution",95,0,95,60,0,6);
 174  fHistchi2ndfvscs->SetXTitle("%cs");
 175  fHistchi2ndfvscs->SetYTitle("Chi2/ndf");
 176
 177
 178 fHistnclus=new TH1D("no of clusters","no of tracks after TPCcluster cut",200,0,200);
 179 fHistchi2ndf= new TH1D("chi_square/ndf","no of tracks after chi_square/ndf cut",10,0,10);
 180
 181  fHistVz = new TH1D("fHistVz","Primary vertex distribution - z coordinate",90,-15,15);
 182  fHistVz->SetXTitle("Vz");
 183  fHistVz->SetYTitle("Entries");
 184
 185
 186   fOutputList->Add(fHistMultV0A);
 187   fOutputList->Add(fHistMultV0C);
 188   fOutputList->Add(fHistMultV0total);
 189   fOutputList->Add(fHistnclus);
 190   fOutputList->Add(fHistchi2ndf);
 191   fOutputList->Add(fHistNchPt);
 192   fOutputList->Add(fHistNchEta);
 193   fOutputList->Add(fHistNchEtaCent);
 194
 195   fOutputList->Add(fHistNchPhi);
 196   fOutputList->Add(fHistVz);
 197   fOutputList->Add(fHistchi2ndfvscs);
 198
 199
 200
 201 TString hname;
 202 TString htitle;
 203
 204   for(Int_t i = 0; i < 25; i++)  {
 205     hname  = "fMult"; hname+=i;
 206     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 207     fMult[i] = new TH1F(hname.Data(),htitle.Data(),5000,0.5,5000.5);
 208     fOutputList->Add(fMult[i]);
 209   }
 210
 211   for(Int_t i = 25; i < 50; i++) {
 212    hname  = "fMult"; hname+=i;
 213     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 214     fMult[i] = new TH1F(hname.Data(),htitle.Data(),3000,0.5,3000.5);
 215     fOutputList->Add(fMult[i]);
 216   }
 217
 218   for(Int_t i = 50; i < 91; i++) {
 219     hname  = "fMult"; hname+=i;
 220     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 221     fMult[i] = new TH1F(hname.Data(),htitle.Data(),2000,0.5,2000.5);
 222     fOutputList->Add(fMult[i]);
 223   }
 224
 225 //****************for 2% centrality bin**************************//
 226 for(Int_t i = 0; i < 12; i++)  {
 227     hname  = "fMultTwo"; hname+=i;
 228     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 229     fMultTwo[i] = new TH1F(hname.Data(),htitle.Data(),5000,0.5,5000.5);
 230     fOutputList->Add(fMultTwo[i]);
 231   }
 232
 233   for(Int_t i = 12; i < 25; i++) {
 234    hname  = "fMultTwo"; hname+=i;
 235     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 236     fMultTwo[i] = new TH1F(hname.Data(),htitle.Data(),3000,0.5,3000.5);
 237     fOutputList->Add(fMultTwo[i]);
 238   }
 239
 240   for(Int_t i = 25; i <= 45; i++) {
 241     hname  = "fMultTwo"; hname+=i;
 242     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 243     fMultTwo[i] = new TH1F(hname.Data(),htitle.Data(),2000,0.5,2000.5);
 244     fOutputList->Add(fMultTwo[i]);
 245   }
 246
 247
 248
 249 //*****************************for 5 % centrality bin********************//
 250 for(Int_t i = 0; i < 5; i++)  {
 251     hname  = "fMultFive"; hname+=i;
 252     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 253     fMultFive[i] = new TH1F(hname.Data(),htitle.Data(),5000,0.5,5000.5);
 254     fOutputList->Add(fMultFive[i]);
 255   }
 256
 257   for(Int_t i = 5; i < 10; i++) {
 258    hname  = "fMultFive"; hname+=i;
 259     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 260     fMultFive[i] = new TH1F(hname.Data(),htitle.Data(),3000,0.5,3000.5);
 261     fOutputList->Add(fMultFive[i]);
 262   }
 263
 264   for(Int_t i = 10; i <= 14; i++) {
 265     hname  = "fMultFive"; hname+=i;
 266     htitle = "Multiplicity in Cent Bin "; htitle+=i;
 267     fMultFive[i] = new TH1F(hname.Data(),htitle.Data(),2000,0.5,2000.5);
 268     fOutputList->Add(fMultFive[i]);
 269   }
 270
 271
 272         }
 273
 274 //________________________________________________________________________
 275 void AliEbyEMultFluctuationTask::UserExec(Option_t *)
 276 {
 277
 278
 279
 280 // Post output data.
 281   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 282   if (!fAOD) {
 283     printf("ERROR: fAOD not available\n");
 284     return;
 285
 286   }
 287
 288   //  Int_t Counter = 0;
 289 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 290 if(!isSelected) return;
 291  fEventCounter->Fill(1);
 292  //Counter++;
 293
 294
 295 // Get the event vertex.
 296         fAODVertex = fAOD->GetPrimaryVertex();
 297
 298         if(TMath::Abs((fAODVertex->GetZ())-((fAOD->GetPrimaryVertexSPD())->GetZ())) > 0.5) return;
 299
 300         if (!fAODVertex) {
 301
 302                 return;
 303         }
 304
 305
 306     if (!SelectEvent(fAODVertex)) return;
 307
 308     //***************VZERO-AMPLITUDE*************************//
 309
 310         AliAODVZERO *aodV0 = fAOD->GetVZEROData();
 311         Float_t fV0Amult = aodV0->GetMTotV0A();
 312         Float_t fV0Cmult = aodV0->GetMTotV0C();
 313         Float_t fV0total = fV0Amult + fV0Cmult;
 314
 315         //*****************SPD-CLUSTERS*************************//
 316
 317
 318
 319                 AliAODHeader *fHeader = fAOD->GetHeader();
 320                 Int_t spdmult0 = fHeader->GetNumberOfITSClusters(0);
 321                 Int_t spdmult1 = fHeader->GetNumberOfITSClusters(1);
 322                 Int_t run = fHeader->GetRunNumber();
 323
 324                 AliAODTracklets *fTracklets = fAOD->GetTracklets();
 325                 Int_t tracklets = fTracklets->GetNumberOfTracklets();
 326
 327
 328                 //********************CENTRALITY******************************************//
 329
 330     AliCentrality *centrality= fAOD->GetCentrality();
 331     Double_t cent=centrality->GetCentralityPercentile("V0M");
 332     //    cout<<"cent= "<<cent<<endl;
 333
 334  if(cent < 0 || cent >= 91) return;
 335
 336 Int_t nCentrality = (Int_t)cent;
 337
 338 Int_t nBin2= (Int_t)( nCentrality/2.0);
 339
 340 Int_t nBin5 = (Int_t)( nCentrality/5.0);
 341  if(nBin5==0){
 342
 343 histcounter->Fill(1);
 344  }
 345
 346
 347
 348 fCentralityCounter->Fill(nCentrality);
 349
 350 fEventCounter->Fill(7);
 351
 352
 353 Double_t Mult=0.0;
 354         Double_t Mult1=0.0;
 355         Double_t Mult2=0.0;
 356
 357         //printf("There are %d tracks in this event\n", fAOD->GetNumberOfTracks());
 358
 359
 360   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
 361   AliAODTrack*  track = fAOD->GetTrack(iTracks);
 362     if (!track) {
 363       printf("ERROR: Could not receive track %d\n", iTracks);
 364       continue;
 365     }
 366
 367
 368
 369
 370 // Find the track classification.
 371 Int_t tracktype = SelectTrack(track);
 372
 373
 374  if (tracktype==0) {
 375    continue;
 376                 }
 377
 378  if(track->GetType()== AliAODTrack::kPrimary)
 379 {
 380   Short_t Charge = track->Charge();
 381         Double_t eta = track->Eta();
 382  if(TMath :: Abs(Charge)==1 && nBin5==0){
 383                        fHistNchEtaCent->Fill(eta);
 384
 385                      }
 386
 387         Double_t pt = track->Pt();
 388         if( pt < 0.3 || pt > 2.0 ) continue;
 389   Double_t nclus = track->GetTPCClusterInfo(2,1);
 390   Double_t chi2ndf = track->Chi2perNDF();
 391   fHistnclus->Fill(nclus);
 392   fHistchi2ndf->Fill(chi2ndf);
 393   fHistchi2ndfvscs->Fill(cent,chi2ndf);
 394
 395                 //Count the charge particles-----
 396                   if(TMath :: Abs(Charge)==1) Mult +=1;
 397                   if(Charge==1) Mult1 +=1;
 398                         if(Charge==-1) Mult2 +=1;
 399
 400
 401         Double_t phi = track->Phi();
 402                      if(TMath :: Abs(Charge)==1) fHistNchPt->Fill(pt);
 403
 404
 405                   if(TMath :: Abs(Charge)==1) fHistNchEta->Fill(eta);
 406                    if(TMath :: Abs(Charge)==1) fHistNchPhi->Fill(phi);
 407
 408
 409  }
 410
 411
 412   } //track loop
 413
 414
 415   //****************Filling Tree*************************//
 416
 417   My_ntuple->Fill(Mult,Mult1,Mult2,nCentrality,fV0total,nBin2,nBin5,spdmult0,spdmult1,tracklets,run);
 418
 419
 420 fMult[nCentrality]->Fill(Mult);
 421
 422 if(nBin2<=45) {
 423 fMultTwo[nBin2]->Fill(Mult);
 424  }
 425 if(nBin5<=14){
 426 fMultFive[nBin5]->Fill(Mult);
 427  }
 428
 429  fHistMultV0A->Fill(fV0Amult);
 430  fHistMultV0C->Fill(fV0Cmult);
 431  fHistMultV0total->Fill(fV0total);
 432
 433   PostData(1, fOutputList);
 434 }
 435 //________________________________________________________________________
 436 void AliEbyEMultFluctuationTask::Terminate(Option_t *)
 437 {
 438
 439
 440   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 441   if (!fOutputList) {
 442     printf("ERROR: Output list not available\n");
 443     return;
 444   }
 445
 446 }