PWG4/JetTasks/AliAnalysisTaskHardSoft.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TH3F.h"
   6 #include "TCanvas.h"
   7 #include "TList.h"
   8 #include "TParticle.h"
   9 #include "TParticlePDG.h"
  10 #include "TProfile.h"
  11 #include "TNtuple.h"
  12 #include "TFile.h"
  13 #include "TRandom.h"
  14
  15 #include "AliAnalysisTask.h"
  16 #include "AliAnalysisManager.h"
  17
  18 #include "AliESDEvent.h"
  19 #include "AliStack.h"
  20 #include "AliMCParticle.h"
  21 #include "AliMCEvent.h"
  22
  23 #include "AliLog.h"
  24 #include "AliESDVertex.h"
  25 #include "AliESDInputHandler.h"
  26 #include "AliESDtrackCuts.h"
  27 #include "AliMultiplicity.h"
  28
  29 #include "AliAnalysisTaskHardSoft.h"
  30 #include "AliExternalTrackParam.h"
  31 #include "AliTrackReference.h"
  32 #include "AliHeader.h"
  33 #include "AliGenEventHeader.h"
  34 #include "AliGenDPMjetEventHeader.h"
  35
  36 // Analysis Task for Hard and Soft event characteristics
  37
  38 // Author: E.Sicking
  39
  40 ClassImp(AliAnalysisTaskHardSoft)
  41
  42 //________________________________________________________________________
  43   AliAnalysisTaskHardSoft::AliAnalysisTaskHardSoft(const char *name)
  44     : AliAnalysisTaskSE(name)
  45     ,fUseMc(false)
  46     ,fUseNeutral(false)
  47     ,fRadiusCut(0.7)
  48     ,fTriggerPtCut(0.7)
  49     ,fAssociatePtCut(0.4)
  50     ,fMap(0x0)
  51     ,fMapLeading(0x0)
  52     ,fCuts(0)
  53     ,fFieldOn(kTRUE)
  54     ,fHists(0)
  55
  56 {
  57   for(Int_t i = 0;i< 2;i++){
  58     fPt[i]                  = 0;
  59     fEta[i]                 = 0;
  60     fPhi[i]                 = 0;
  61     fEtaPhi[i]              = 0;
  62     fEtaPhiLeading[i]       = 0;
  63     fNch[i]                 = 0;
  64     fPtLeading[i]           = 0;
  65     fPtLeadingNch[i]        = 0;
  66     fPtSumNch[i]            = 0;
  67     fPtAvNch[i]             = 0;
  68     fDPhiLeading[i]         = 0;
  69     fRadiusLeading[i]       = 0;
  70     fDPhiLeadingR[i]        = 0;
  71     fRadiusLeadingR[i]      = 0;
  72     fDPhiLeadingRS[i]       = 0;
  73     fRadiusLeadingRS[i]     = 0;
  74     fNchAssInR[i]           = 0;
  75     fTrigger[i]             = 0;
  76
  77     for(Int_t j=0;j<100;j++){
  78       fDPhiLeadingNchBin[i][j]   = 0;
  79     }
  80
  81     for(Int_t j=0;j<2;j++){
  82       fNchHardSoft[i][j]   = 0;
  83       fPtHardSoft[i][j]   = 0;
  84       fPtAvHardSoftNch[i][j]   = 0;
  85     }
  86   }
  87   DefineOutput(1,  TList::Class());
  88 }
  89
  90
  91 //________________________________________________________________________
  92 void AliAnalysisTaskHardSoft::UserCreateOutputObjects()
  93 {
  94   // Create histograms
  95   // Called once
  96
  97   Bool_t oldStatus = TH1::AddDirectoryStatus();
  98   TH1::AddDirectory(kFALSE);
  99
 100   fHists = new TList();
 101
 102
 103
 104   TString labels[2];
 105   labels[0]="MC";
 106   labels[1]="ESD";
 107
 108   for(Int_t i=0;i<2;i++){
 109     fPt[i]                 = new TH1F(Form("fPt%s",labels[i].Data()),
 110                                       Form("fPt%s",labels[i].Data()) ,
 111                                       500, 0., 50);
 112     fEta[i]                = new TH1F (Form("fEta%s",labels[i].Data()),
 113                                        Form("fEta%s",labels[i].Data()) ,
 114                                        100, -1., 1);
 115     fPhi[i]                = new TH1F(Form("fPhi%s",labels[i].Data()),
 116                                       Form("fPhi%s",labels[i].Data()) ,
 117                                       360, 0.,2*TMath::Pi());
 118     fEtaPhi[i]             = new TH2F(Form("fEtaPhi%s",labels[i].Data()),
 119                                       Form("fEtaPhi%s",labels[i].Data()) ,
 120                                       100,-1.,1.,
 121                                       360, 0.,2*TMath::Pi());
 122     fEtaPhiLeading[i]      = new TH2F(Form("fEtaPhiLeading%s",labels[i].Data()),
 123                                       Form("fEtaPhiLeading%s",labels[i].Data()) ,
 124                                       100,-1.,1.,
 125                                       360, 0.,2*TMath::Pi());
 126     fNch[i]                = new TH1F(Form("fNch%s",labels[i].Data()),
 127                                       Form("fNch%s",labels[i].Data()) ,
 128                                       250, -0.5, 249.5);
 129     fPtLeading[i]          = new TH1F(Form("fPtLeading%s",labels[i].Data()),
 130                                       Form("fPtLeading%s",labels[i].Data()) ,
 131                                       500, 0., 50);
 132     fPtLeadingNch[i]       = new TProfile(Form("fPtLeadingNch%s",labels[i].Data()),
 133                                           Form("fPtLeadingNch%s",labels[i].Data()) ,
 134                                           250, -0.5, 249.5);
 135     fPtSumNch[i]           = new TProfile(Form("fPtSumNch%s",labels[i].Data()),
 136                                           Form("fPtSumNch%s",labels[i].Data()) ,
 137                                           250, -0.5, 249.5);
 138     fPtAvNch[i]            = new TProfile(Form("fPtAvNch%s",labels[i].Data()),
 139                                           Form("fPtAvNch%s",labels[i].Data()) ,
 140                                           250, -0.5, 249.5);
 141     fDPhiLeading[i]        = new TH1F(Form("fDPhiLeading%s",labels[i].Data()),
 142                                       Form("fDPhiLeading%s",labels[i].Data()) ,
 143                                       180, 0., TMath::Pi());
 144     fRadiusLeading[i]      = new TH1F(Form("fRadiusLeading%s",labels[i].Data()),
 145                                       Form("fRadiusLeading%s",labels[i].Data()) ,
 146                                       180, 0., 2*TMath::Pi());
 147     fDPhiLeadingR[i]       = new TH1F(Form("fDPhiLeadingR%s",labels[i].Data()),
 148                                       Form("fDPhiLeadingR%s",labels[i].Data()) ,
 149                                       180, 0., TMath::Pi());
 150     fRadiusLeadingR[i]     = new TH1F(Form("fRadiusLeadingR%s",labels[i].Data()),
 151                                       Form("fRadiusLeadingR%s",labels[i].Data()) ,
 152                                       180, 0., 2*TMath::Pi());
 153     fDPhiLeadingRS[i]      = new TH1F(Form("fDPhiLeadingRS%s",labels[i].Data()),
 154                                       Form("fDPhiLeadingRS%s",labels[i].Data()) ,
 155                                       180, 0., TMath::Pi());
 156     fRadiusLeadingRS[i]    = new TH1F(Form("fRadiusLeadingRS%s",labels[i].Data()),
 157                                       Form("fRadiusLeadingRS%s",labels[i].Data()) ,
 158                                       180, 0., 2*TMath::Pi());
 159     fNchAssInR[i]          = new TProfile(Form("fNchAssInR%s",labels[i].Data()),
 160                                           Form("fNchAssInR%s",labels[i].Data()) ,
 161                                           250, -0.5, 249.5);
 162     fTrigger[i]            = new TH1F(Form("fTrigger%s",labels[i].Data()),
 163                                       Form("fTrigger%s",labels[i].Data()) ,
 164                                       250, -0.5, 249.5);
 165     for(Int_t j=0;j<100;j++){
 166       fDPhiLeadingNchBin[i][j]     = new TH1F(Form("fDPhiLeadingNchBin%s%02d",labels[i].Data(),j),
 167                                               Form("fDPhiLeadingNchBin%s%02d",labels[i].Data(),j) ,
 168                                               180, 0., TMath::Pi());
 169     }
 170
 171     for(Int_t j=0;j<2;j++){
 172       fNchHardSoft[i][j]     = new TH1F(Form("fNchHardSoft%s%d",labels[i].Data(),j),
 173                                         Form("fNchHardSoft%s%d",labels[i].Data(),j) ,
 174                                         250, -0.5, 249.5);
 175       fPtHardSoft[i][j]     = new TH1F(Form("fPtHardSoft%s%d",labels[i].Data(),j),
 176                                        Form("fPtHardSoft%s%d",labels[i].Data(),j) ,
 177                                        500, 0., 50.0);
 178       fPtAvHardSoftNch[i][j]     = new TProfile(Form("fPtAvHardSoftNch%s%d",labels[i].Data(),j),
 179                                                 Form("fPtAvHardSoftNch%s%d",labels[i].Data(),j) ,
 180                                                 250, -0.5, 249.5);
 181     }
 182
 183   }
 184
 185
 186   fHists->SetOwner();
 187   for(Int_t i=0;i<2;i++){
 188     fHists->Add(fPt[i]);
 189     fHists->Add(fEta[i]);
 190     fHists->Add(fPhi[i]);
 191     fHists->Add(fEtaPhi[i]);
 192     fHists->Add(fEtaPhiLeading[i]);
 193     fHists->Add(fNch[i]);
 194     fHists->Add(fPtLeading[i]);
 195     fHists->Add(fPtLeadingNch[i]);
 196     fHists->Add(fPtSumNch[i]);
 197     fHists->Add(fPtAvNch[i]);
 198     fHists->Add(fDPhiLeading[i]);
 199     fHists->Add(fRadiusLeading[i]);
 200     fHists->Add(fDPhiLeadingR[i]);
 201     fHists->Add(fRadiusLeadingR[i]);
 202     fHists->Add(fDPhiLeadingRS[i]);
 203     fHists->Add(fRadiusLeadingRS[i]);
 204     fHists->Add(fNchAssInR[i]);
 205     fHists->Add(fTrigger[i]);
 206
 207     for(Int_t j=0;j<100;j++){
 208       fHists->Add(fDPhiLeadingNchBin[i][j]);
 209     }
 210
 211     for(Int_t j=0;j<2;j++){
 212       fHists->Add(fNchHardSoft[i][j]);
 213       fHists->Add(fPtHardSoft[i][j]);
 214       fHists->Add(fPtAvHardSoftNch[i][j]);
 215     }
 216   }
 217
 218
 219   TH1::AddDirectory(oldStatus);
 220 }
 221
 222 //__________________________________________________________
 223
 224 void AliAnalysisTaskHardSoft::UserExec(Option_t *)
 225 {
 226   Int_t nentries[2]         = {0,0};  // tracks in event (before selection)
 227   Int_t nTracksAll[2]       = {0,0};  // accepted tracks in event
 228   Int_t nTracksAssociate[2] = {0,0};  // accepted tracks in event within R=fRadiusCut around leading track
 229
 230   Float_t pt[2]    = {0.,0.}; // pt
 231   Float_t eta[2]   = {0.,0.}; // eta
 232   Float_t phi[2]   = {0.,0.}; // phi
 233   Float_t ptSum[2] = {0.,0.}; // pt sum
 234   Float_t ptAv[2]  = {0.,0.}; // average pt
 235
 236
 237   Float_t ptLeading[2]     = {0.,0.}; // pt leading
 238   Float_t etaLeading[2]    = {0.,0.}; // eta leading
 239   Float_t phiLeading[2]    = {0.,0.}; // phi leading
 240
 241   Float_t ptOthers[2]     = {0.,0.}; // pt others // for second track loop around leading track
 242   Float_t etaOthers[2]    = {0.,0.}; // eta others
 243   Float_t phiOthers[2]    = {0.,0.}; // phi others
 244
 245   Double_t etaLeadingRandom[2]   = {0.,0.}; // eta leading for random particle position (from fMapLeading)
 246   Double_t phiLeadingRandom[2]   = {0.,0.}; // phi leading  "
 247   Double_t etaOthersRandom[2]    = {0.,0.}; // eta others  for random particle position (from fMap)
 248   Double_t phiOthersRandom[2]    = {0.,0.}; // phi others   "
 249
 250   Int_t fEventType[2]={1,1}; //hard=0, soft=1 -> set to hard if trigger and associate particle
 251                              //are within R=0.7
 252
 253   //get event and vertex cut :(MC and ESD)
 254   //---------------------
 255   //MC
 256   //---------------------
 257   AliStack *stack = 0x0; // needed for MC studies
 258   Float_t vzMC=0.;       // MC vertex position in z
 259   if(fUseMc==true){
 260     stack = MCEvent()->Stack();
 261     AliGenEventHeader*  header = MCEvent()->GenEventHeader();
 262     TArrayF mcV;
 263     header->PrimaryVertex(mcV);
 264     vzMC = mcV[2];
 265   }
 266   //---------------------
 267   //ESD
 268   //---------------------
 269   AliVEvent* event = InputEvent();
 270   if (!event) {
 271     Printf("ERROR: Could not retrieve event");
 272     return;
 273   }
 274   if(Entry()==0){
 275     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 276     if(esd)Printf("We are reading from ESD");
 277   }
 278   const AliVVertex* vertex = event->GetPrimaryVertex();
 279   Float_t vz = vertex->GetZ();
 280   //---------------------
 281
 282
 283
 284
 285   //Selection of particle(mcesd=0) or esdtracks(mcesd=1)
 286   for(Int_t mcesd=0;mcesd<2;mcesd++){
 287     if(mcesd==0){
 288       if(fUseMc==false) //MC part can be switched off for real data by function SetUserMc(kFALSE)
 289         continue;
 290     }
 291
 292     // vertex cut and number of particles/tracks per event
 293     //---------------------------------------
 294     if(mcesd==0){//mc particles
 295       if (TMath::Abs(vzMC) > 10.) return;
 296       nentries[mcesd]=MCEvent()->GetNumberOfTracks();
 297     }
 298     else{// esd tracks
 299       if (TMath::Abs(vz) > 10.) return;
 300       nentries[mcesd]=event->GetNumberOfTracks();
 301     }//---------------------------------------
 302
 303
 304
 305
 306     // arrays for leading track determination (done with TMath::Sort of Array)
 307     Float_t * ptArray = new Float_t[nentries[mcesd]];
 308
 309     //array of track properties
 310     Float_t * etaArray = new Float_t[nentries[mcesd]];
 311     Float_t * phiArray = new Float_t[nentries[mcesd]];
 312
 313     Int_t *pindex  = new Int_t[nentries[mcesd]];
 314     for (Int_t i = 0; i < nentries[mcesd]; i++) {
 315       ptArray[i]=0.;
 316       etaArray[i]=0.;
 317       phiArray[i]=0.;
 318       pindex[i]=0;
 319     }
 320
 321
 322
 323
 324     //first track loop
 325     for (Int_t iTrack = 0; iTrack < nentries[mcesd]; iTrack++) {
 326
 327       pt[mcesd]  = 0.;
 328       eta[mcesd] = 0.;
 329       phi[mcesd] = 0;
 330
 331       //get properties of mc particle
 332       if(mcesd==0){//mc
 333         AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
 334         // Primaries only
 335         if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
 336         if(!fUseNeutral)if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
 337         //same cuts as on ESDtracks
 338         if(TMath::Abs(mcP->Eta())>0.9)continue;
 339         if(mcP->Pt()<0.2)continue;
 340         if(mcP->Pt()>200)continue;
 341
 342         pt[mcesd]  = mcP->Pt();
 343         eta[mcesd] = mcP->Eta();
 344         phi[mcesd] = mcP->Phi();
 345       }
 346
 347       //get properties of esdtracks
 348       else{//esd
 349         AliVParticle *track = event->GetTrack(iTrack);
 350         if (!track) {
 351           Printf("ERROR: Could not receive track %d", iTrack);
 352           continue;
 353         }
 354         AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);
 355         if(!esdtrack)continue;
 356         if (!fCuts->AcceptTrack(esdtrack)) continue;
 357
 358         pt[mcesd]=esdtrack->Pt();
 359         eta[mcesd]=esdtrack->Eta();
 360         phi[mcesd]=esdtrack->Phi();
 361       }
 362
 363       ptArray[nTracksAll[mcesd]]    = pt[mcesd];  // fill pt array
 364       etaArray[nTracksAll[mcesd]]   = eta[mcesd]; // fill eta array
 365       phiArray[nTracksAll[mcesd]++] = phi[mcesd]; // count tracks and fill phi array
 366
 367       fPt[mcesd]     -> Fill(pt[mcesd]);
 368       fEta[mcesd]    -> Fill(eta[mcesd]);
 369       fPhi[mcesd]    -> Fill(phi[mcesd]);
 370       fEtaPhi[mcesd] -> Fill(eta[mcesd],phi[mcesd]);
 371
 372
 373     }//end first track loop
 374
 375     fNch[mcesd]    -> Fill(nTracksAll[mcesd]);
 376
 377
 378     //find leading pt tracks
 379     if(nentries[mcesd]>0) TMath::Sort(nentries[mcesd], ptArray, pindex, kTRUE);
 380     //for(Int_t i=0;i<nTracksAll[mcesd];i++){//show just the filled entries, skip empty ones.
 381     //     printf("%d:  pt = %f, number %i \n",mcesd, ptArray[pindex[i]],i);
 382     //}
 383
 384
 385     if(nTracksAll[mcesd]>0){
 386       fPtLeading[mcesd]->Fill(ptArray[pindex[0]]);  //first entry in array is highest
 387       fPtLeadingNch[mcesd]->Fill(nTracksAll[mcesd],ptArray[pindex[0]]);
 388
 389       for(Int_t i=0;i<nTracksAll[mcesd];i++){      //calculate ptsum
 390         ptSum[mcesd]+=ptArray[pindex[i]];
 391       }
 392       ptAv[mcesd]=ptSum[mcesd]/nTracksAll[mcesd];  //calculate <pt>
 393
 394       fPtSumNch[mcesd]->Fill(nTracksAll[mcesd],ptSum[mcesd]);
 395       fPtAvNch[mcesd]->Fill(nTracksAll[mcesd],ptAv[mcesd]);
 396     }
 397
 398
 399     if(nTracksAll[mcesd]>1){ // require at least two tracks (leading and prob. accosicates)
 400
 401       //Leading track properties
 402       ptLeading[mcesd]  = ptArray[pindex[0]];
 403       etaLeading[mcesd] = etaArray[pindex[0]];
 404       phiLeading[mcesd] = phiArray[pindex[0]];
 405
 406       fEtaPhiLeading[mcesd] -> Fill(etaLeading[mcesd],phiLeading[mcesd]);
 407
 408       fMapLeading->GetRandom2(etaLeadingRandom[mcesd],phiLeadingRandom[mcesd]);
 409
 410       //second track loop for event propoerties around leading tracks with pt>triggerPtCut
 411       //loop only over already accepted tracks except leading track
 412       if(ptLeading[mcesd]>fTriggerPtCut){
 413
 414         fTrigger[mcesd]->Fill(nTracksAll[mcesd]); // how often is there a trigger particle at a certain Nch bin
 415
 416         for (Int_t iTrack = 1; iTrack < nTracksAll[mcesd]; iTrack++) { // starting at second highest track
 417
 418           ptOthers[mcesd]   = ptArray[pindex[iTrack]];
 419           etaOthers[mcesd]  = etaArray[pindex[iTrack]];
 420           phiOthers[mcesd]  = phiArray[pindex[iTrack]];
 421
 422           fMap->GetRandom2(etaOthersRandom[mcesd],phiOthersRandom[mcesd]);
 423
 424           if(ptOthers[mcesd]>fAssociatePtCut){ // only tracks which fullfill associate pt cut
 425
 426             //1. real data
 427
 428             Float_t dPhi=TMath::Abs(phiOthers[mcesd]-phiLeading[mcesd]);
 429             if(dPhi>TMath::Pi())      dPhi=2*TMath::Pi()-dPhi;
 430             Float_t dEta=etaOthers[mcesd]-etaLeading[mcesd];
 431
 432             Float_t radius=TMath::Sqrt(dPhi*dPhi+dEta*dEta);
 433             fRadiusLeading[mcesd]->Fill(radius);
 434             fDPhiLeading[mcesd]->Fill(dPhi);
 435             if(nTracksAll[mcesd]<100)fDPhiLeadingNchBin[mcesd][nTracksAll[mcesd]]->Fill(dPhi);
 436
 437             if(radius<fRadiusCut){
 438               fEventType[mcesd]=0;
 439               nTracksAssociate[mcesd]++;
 440             }
 441
 442             //2. random position
 443             Float_t dPhiR=TMath::Abs(phiOthersRandom[mcesd]-phiLeadingRandom[mcesd]);
 444             if(dPhiR>TMath::Pi())      dPhiR=dPhiR-2*TMath::Pi();
 445             Float_t dEtaR=etaOthersRandom[mcesd]-etaLeadingRandom[mcesd];
 446
 447             Float_t radiusR=TMath::Sqrt(dPhiR*dPhiR+dEtaR*dEtaR);
 448             fRadiusLeadingR[mcesd]->Fill(radiusR);
 449             fDPhiLeadingR[mcesd]->Fill(dPhiR);
 450
 451
 452
 453             //3. random position of leading particle
 454             Float_t dPhiRS=TMath::Abs(phiOthers[mcesd]-phiLeadingRandom[mcesd]);
 455             if(dPhiRS>TMath::Pi())      dPhiRS=dPhiRS-2*TMath::Pi();
 456             Float_t dEtaRS=etaOthers[mcesd]-etaLeadingRandom[mcesd];
 457
 458             Float_t radiusRS=TMath::Sqrt(dPhiRS*dPhiRS+dEtaRS*dEtaRS);
 459             fRadiusLeadingRS[mcesd]->Fill(radiusRS);
 460             fDPhiLeadingRS[mcesd]->Fill(dPhiRS);
 461           }
 462
 463
 464
 465         }
 466         //fill histogram with number of tracks (pt>fAssociatePtCut) around leading track
 467         fNchAssInR[mcesd]->Fill(nTracksAll[mcesd],nTracksAssociate[mcesd]);
 468
 469       }
 470     }
 471
 472     fNchHardSoft[mcesd][fEventType[mcesd]]->Fill(nTracksAll[mcesd]);
 473
 474     for(Int_t i=0;i<nTracksAll[mcesd];i++){
 475       fPtHardSoft[mcesd][fEventType[mcesd]]->Fill(ptArray[i]);
 476     }
 477
 478     fPtAvHardSoftNch[mcesd][fEventType[mcesd]]->Fill(nTracksAll[mcesd],ptAv[mcesd]);
 479
 480
 481   }//double loop over mcP and ESD
 482
 483
 484
 485   // Post output data.
 486   PostData(1, fHists);
 487 }
 488
 489
 490
 491
 492
 493 //________________________________________________________________________
 494 void AliAnalysisTaskHardSoft::Terminate(Option_t *)
 495 {
 496
 497
 498 }
 499
 500
 501
 502
 503