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