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