PWG4/JetTasks/AliAnalysisTaskQGSep.cxx

   1 #include <TChain.h>
   2 #include <TList.h>
   3
   4 #include <TTree.h>
   5 #include <TFile.h>
   6 #include <TH1.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TProfile.h>
  10 #include <TCanvas.h>
  11
  12 #include "AliVEvent.h"
  13 #include "AliVParticle.h"
  14
  15 #include "AliESDEvent.h"
  16 #include "AliESDtrack.h"
  17
  18 #include "AliAODEvent.h"
  19 #include "AliAODTrack.h"
  20
  21 #include "AliMCEvent.h"
  22 #include "AliMCParticle.h"
  23 #include "AliAnalysisManager.h"
  24 #include "AliMCEvent.h"
  25 #include "TParticle.h"
  26 #include "AliStack.h"
  27
  28 #include "AliAODEvent.h"
  29 #include "AliAODVertex.h"
  30 #include "AliAODHandler.h"
  31 #include "AliAODTrack.h"
  32 #include "AliAODJet.h"
  33 #include "AliAODMCParticle.h"
  34
  35 #include "AliAnalysisTaskQGSep.h"
  36 #include "AliAnalysisHelperJetTasks.h"
  37 ClassImp(AliAnalysisTaskQGSep)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskQGSep::AliAnalysisTaskQGSep(const char *name)
  41   : AliAnalysisTaskSE(name),
  42     fBranchRec("jets"),
  43     fUseMC(kFALSE),
  44     fUseAOD(kFALSE),
  45     fXsection(0),
  46     fWeight(0),
  47     fMyAODEvent(0),
  48     fOutputList(0),
  49     fpHistPtAvEQ(0),
  50     fpHistPtAvEG(0),
  51     fpHistDrEQ(0),
  52     fpHistDrEG(0),
  53     fpHistDrE(0),
  54     fpHistPtAvE(0),
  55     fpHistDrE3(0),
  56     fpHistPtAvE3(0)
  57
  58 {
  59   // Constructor
  60
  61   DefineOutput(1, TList::Class());
  62 }
  63
  64 //________________________________________________________________________
  65 void AliAnalysisTaskQGSep::UserCreateOutputObjects()
  66 {
  67   // Create histograms
  68   // Called once
  69
  70   fOutputList = new TList();
  71
  72   if(fUseMC){
  73
  74     //histos for quarks
  75     fpHistPtAvEQ = new TProfile("fpHistPtAvEQ", "", 100, 0, 100, 0, 10);
  76     fOutputList->Add(fpHistPtAvEQ);
  77     fpHistDrEQ = new TProfile("fpHistDrEQ", "", 100, 0, 100, 0, 0.5);
  78     fOutputList->Add(fpHistDrEQ);
  79
  80     //histos for gluons
  81     fpHistPtAvEG = new TProfile("fpHistPtAvEG", "", 100, 0, 100, 0., 10.);
  82     fOutputList->Add(fpHistPtAvEG);
  83     fpHistDrEG = new TProfile("fpHistDrEG", "", 100, 0, 100, 0, 0.5);
  84     fOutputList->Add(fpHistDrEG);
  85   }
  86
  87   //histos for full spetra, no separation
  88   fpHistDrE = new TProfile("fpHistDrE", "", 100, 0, 100, 0, 0.5);
  89   fOutputList->Add(fpHistDrE);
  90   fpHistPtAvE = new TProfile("fpHistPtAvE", "", 100, 0, 100, 0., 10.);
  91   fOutputList->Add(fpHistPtAvE);
  92   fpHistDrE3 = new TProfile("fpHistDrE3", "", 100, 0, 100, 0, 0.5);
  93   fOutputList->Add(fpHistDrE3);
  94   fpHistPtAvE3 = new TProfile("fpHistPtAvE3", "", 100, 0, 100, 0., 10.);
  95   fOutputList->Add(fpHistPtAvE3);
  96
  97   for(Int_t i = 0; i < fOutputList->GetEntries(); i++){
  98     TH1 * h = dynamic_cast<TH1*>(fOutputList->At(i));
  99     if(h) h->Sumw2();
 100   }
 101
 102   if(fDebug)Printf("~# QGSep User objects created");
 103 }
 104
 105 //__________________________________________________________________
 106 Bool_t AliAnalysisTaskQGSep::Notify()
 107 {
 108
 109   if(fUseMC){
 110     //
 111     // Implemented Notify() to read the cross sections
 112     // and number of trials from pyxsec.root
 113     //
 114
 115
 116     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 117     Float_t xsection = 0;
 118     Float_t ftrials  = 1;
 119
 120     Float_t fAvgTrials = 1;
 121     if(tree){
 122       TFile *curfile = tree->GetCurrentFile();
 123       if (!curfile) {
 124         Error("Notify","No current file");
 125         return kFALSE;
 126       }
 127       AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 128       // construct a poor man average trials
 129       Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 130       if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; // CKB take this into account for normalisation
 131     }
 132
 133     if(xsection>0){
 134       fXsection  = xsection;
 135       fWeight = fXsection/fAvgTrials;
 136     }
 137     else fWeight = 1;
 138     return kTRUE;
 139   }
 140
 141   else{
 142     fWeight = 1;
 143     return kTRUE;
 144   }
 145
 146 }
 147
 148 //________________________________________________________________________
 149 void AliAnalysisTaskQGSep::UserExec(Option_t *)
 150 {
 151   // Main loop
 152   // Called for each event
 153
 154   if(fUseAOD){
 155     AliVEvent *event = InputEvent();
 156     if (!event) {
 157       Error("UserExec", "Could not retrieve event");
 158       return;
 159     }
 160
 161   fMyAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
 162   }
 163
 164   else{
 165         //  assume that the AOD is in the general output...
 166     fMyAODEvent  = AODEvent();
 167     if(!fMyAODEvent){
 168       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 169       return;
 170     }
 171   }
 172
 173   if (fMyAODEvent){
 174     if(fUseMC)
 175       LoopAODMC();
 176     else
 177       LoopAOD();
 178   }
 179
 180   // Post output data.
 181   PostData(1, fOutputList);
 182 }
 183
 184 //________________________________________________________________________
 185 void AliAnalysisTaskQGSep::Terminate(Option_t*)
 186 {
 187   // Draw result to the screen
 188   // Called once at the end of the query
 189
 190   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 191   if (!fOutputList) {
 192     Error("Terminate","fOutputList not available");
 193     return;
 194   }
 195
 196 }
 197
 198 //________________________________________________________________________
 199 void AliAnalysisTaskQGSep::LoopAOD(){
 200   AliAODJet recJets[4];
 201   AliAODJet jets[4];
 202   AliAODJet rJets[4];
 203   Int_t nRecJets = 0;
 204
 205   //array of reconstructed jets from the AOD input
 206   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fMyAODEvent->FindListObject(fBranchRec.Data()));
 207   if(!aodRecJets){
 208     PostData(1, fOutputList);
 209     return;
 210   }
 211
 212   // reconstructed jets
 213   nRecJets = aodRecJets->GetEntries();
 214   if(fDebug)Printf("--- Jets found in bRec: %d", nRecJets);
 215   nRecJets = TMath::Min(nRecJets, 4);
 216   for(int ir = 0;ir < nRecJets;++ir)
 217     {
 218       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 219       if(!tmp)continue;
 220       jets[ir] = *tmp;
 221     }
 222
 223   Int_t counter = 0;
 224   Int_t tag = 0;
 225   Int_t nGenSel = 0;
 226
 227   TLorentzVector v[4];
 228   Double_t eSum = 0.;
 229   Double_t pxSum = 0.;
 230   Double_t pySum = 0.;
 231   Double_t pzSum = 0.;
 232
 233   for(Int_t i = 0; i < nRecJets; i++)
 234     {
 235       if(nRecJets == 1)
 236         {
 237           rJets[nGenSel] = jets[i];
 238           v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
 239           eSum += jets[i].E();
 240           pxSum += jets[i].Px();
 241           pySum += jets[i].Py();
 242           pzSum += jets[i].Pz();
 243           nGenSel++;
 244         }
 245       else
 246         {
 247           counter = 0;
 248           tag = 0;
 249           for(Int_t j = 0; j < nRecJets; j++)
 250             {
 251               if(i!=j)
 252                 {
 253                   Double_t dRij = jets[i].DeltaR(&jets[j]);
 254                   counter++;
 255                   if(dRij > 2*0.4) tag++;
 256                 }
 257             }
 258           if(counter!=0)
 259             {
 260               if(tag/counter == 1)
 261                 {
 262                   rJets[nGenSel] = jets[i];
 263                   v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
 264                   eSum += jets[i].E();
 265                   pxSum += jets[i].Px();
 266                   pySum += jets[i].Py();
 267                   pzSum += jets[i].Pz();
 268                   nGenSel++;
 269                 }
 270             }
 271         }
 272     }
 273
 274   nRecJets = nGenSel;
 275
 276   if(nRecJets == 0){
 277     PostData(1, fOutputList);
 278     return;
 279   }
 280
 281   TLorentzVector vB;
 282   vB.SetPxPyPzE(pxSum, pySum, pzSum, eSum);
 283   Double_t e[4];
 284   for(Int_t i = 0; i < nRecJets; i++){
 285     v[i].Boost(-vB.Px()/vB.E(),-vB.Py()/vB.E(),-vB.Pz()/vB.E());
 286     e[i] = v[i].E();
 287   }
 288
 289   Int_t idxj[4];
 290   TMath::Sort(nRecJets, e, idxj);
 291   for(Int_t i = 0; i < nRecJets; i++){
 292     recJets[i] = rJets[idxj[i]];
 293   }
 294
 295
 296   for(Int_t iJ = 0; iJ < nRecJets; iJ++){
 297     Double_t pTsum = 0.;
 298     TRefArray * tra = dynamic_cast<TRefArray*>(recJets[iJ].GetRefTracks());
 299     if(!tra) continue;
 300     Int_t nAODtracks = TMath::Min(1000, tra->GetEntries());
 301     Double_t dR[1000];
 302     for(Int_t iT = 0; iT < nAODtracks; iT++){
 303       AliAODTrack * jetTrack = dynamic_cast<AliAODTrack*>(tra->At(iT));
 304       if(!jetTrack) continue;
 305       pTsum += jetTrack->Pt();
 306       dR[iT] = recJets[iJ].DeltaR(jetTrack);
 307     }
 308
 309
 310     fpHistPtAvE->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
 311
 312     if(iJ > 1) //fill mulit-jet histo
 313       fpHistPtAvE3->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
 314
 315     Int_t idxAOD[1000];
 316     TMath::Sort(nAODtracks, dR, idxAOD, kFALSE);
 317
 318     Double_t pTsum90Inv=0.;
 319     for(Int_t iT = 0; iT < nAODtracks; iT++){
 320       AliAODTrack * track = dynamic_cast<AliAODTrack*>(tra->At(idxAOD[iT]));
 321       if(!track) continue;
 322       pTsum90Inv += track->Pt();
 323       if(pTsum90Inv >= 0.9*pTsum){
 324         Double_t deltaR = recJets[iJ].DeltaR(track);
 325
 326         fpHistDrE->Fill(recJets[iJ].E(), deltaR, fWeight);
 327
 328         if(iJ > 1)
 329           fpHistDrE3->Fill(recJets[iJ].E(), deltaR, fWeight);
 330
 331         break;
 332       }
 333     }
 334   }
 335 }
 336
 337 //__________________________________________________________________
 338 void AliAnalysisTaskQGSep::LoopAODMC(){
 339
 340   AliAODJet recJets[4];
 341   AliAODJet jets[4];
 342   AliAODJet rJets[4];
 343   Int_t nRecJets = 0;
 344
 345   //array of reconstructed jets from the AOD input
 346   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fMyAODEvent->FindListObject(fBranchRec.Data()));
 347   if(!aodRecJets){
 348     PostData(1, fOutputList);
 349     return;
 350   }
 351
 352   // reconstructed jets
 353   nRecJets = aodRecJets->GetEntries();
 354   if(fDebug)Printf("--- Jets found in bRec: %d", nRecJets);
 355   nRecJets = TMath::Min(nRecJets, 4);
 356   for(int ir = 0;ir < nRecJets;++ir)
 357     {
 358       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 359       if(!tmp)continue;
 360       jets[ir] = *tmp;
 361     }
 362
 363   Int_t counter = 0;
 364   Int_t tag = 0;
 365   Int_t nGenSel = 0;
 366
 367   TLorentzVector v[4];
 368   Double_t eSum = 0.;
 369   Double_t pxSum = 0.;
 370   Double_t pySum = 0.;
 371   Double_t pzSum = 0.;
 372
 373   for(Int_t i = 0; i < nRecJets; i++)
 374     {
 375       if(nRecJets == 1)
 376         {
 377           rJets[nGenSel] = jets[i];
 378           v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
 379           eSum += jets[i].E();
 380           pxSum += jets[i].Px();
 381           pySum += jets[i].Py();
 382           pzSum += jets[i].Pz();
 383           nGenSel++;
 384         }
 385       else
 386         {
 387           counter = 0;
 388           tag = 0;
 389           for(Int_t j = 0; j < nRecJets; j++)
 390             {
 391               if(i!=j)
 392                 {
 393                   Double_t dRij = jets[i].DeltaR(&jets[j]);
 394                   counter++;
 395                   if(dRij > 2*0.4) tag++;
 396                 }
 397             }
 398           if(counter!=0)
 399             {
 400               if(tag/counter == 1)
 401                 {
 402                   rJets[nGenSel] = jets[i];
 403                   v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
 404                   eSum += jets[i].E();
 405                   pxSum += jets[i].Px();
 406                   pySum += jets[i].Py();
 407                   pzSum += jets[i].Pz();
 408                   nGenSel++;
 409                 }
 410             }
 411         }
 412     }
 413
 414   nRecJets = nGenSel;
 415
 416   if(nRecJets == 0){
 417     PostData(1, fOutputList);
 418     return;
 419    }
 420
 421   TLorentzVector vB;
 422   vB.SetPxPyPzE(pxSum, pySum, pzSum, eSum);
 423   Double_t e[4];
 424   for(Int_t i = 0; i < nRecJets; i++){
 425     v[i].Boost(-vB.Px()/vB.E(),-vB.Py()/vB.E(),-vB.Pz()/vB.E());
 426     e[i] = v[i].E();
 427   }
 428
 429   Int_t idxj[4];
 430   TMath::Sort(nRecJets, e, idxj);
 431   for(Int_t i = 0; i < nRecJets; i++){
 432     recJets[i] = rJets[idxj[i]];
 433   }
 434
 435   Int_t nMCtracks = 0;
 436   TClonesArray *tca = dynamic_cast<TClonesArray*>(fMyAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 437   if(!tca) return;
 438   nMCtracks = TMath::Min(tca->GetEntries(), 10000);
 439
 440   //sort AliAODMCParticles in pT
 441   Double_t pTMC[10000];
 442   for(Int_t iMC = 0; iMC < nMCtracks; iMC++){
 443     AliAODMCParticle *partMC = dynamic_cast<AliAODMCParticle*>(tca->At(iMC));
 444     if(!partMC) continue;
 445     pTMC[iMC]=partMC->Pt();
 446   }
 447
 448   Int_t idxMC[10000];
 449   TMath::Sort(nMCtracks, pTMC, idxMC);
 450
 451
 452   Int_t flagQ[4], flagG[4];
 453   for(Int_t iJ = 0; iJ < nRecJets; iJ++){
 454     //flag jet as q/g
 455
 456     //look for highest momentum parton in the jet cone
 457     for(Int_t iMC = 0; iMC < nMCtracks; iMC++){
 458       AliAODMCParticle *partMC = dynamic_cast<AliAODMCParticle*>(tca->At(idxMC[iMC]));
 459       if(!partMC) continue;
 460       Double_t r = recJets[iJ].DeltaR(partMC);
 461       if(r < 0.4){
 462         if(TMath::Abs(partMC->GetPdgCode()) < 9)
 463           flagQ[iJ] = 1;
 464         if(TMath::Abs(partMC->GetPdgCode()) == 21)
 465           flagG[iJ] = 1;
 466         break;
 467       }
 468     }
 469
 470     Double_t pTsum = 0.;
 471     TRefArray * tra = dynamic_cast<TRefArray*>(recJets[iJ].GetRefTracks());
 472     if(!tra) continue;
 473     Int_t nAODtracks = TMath::Min(1000, tra->GetEntries());
 474     Double_t dR[1000];
 475     for(Int_t iT = 0; iT < nAODtracks; iT++){
 476       AliAODTrack * jetTrack = dynamic_cast<AliAODTrack*>(tra->At(iT));
 477       if(!jetTrack) continue;
 478       pTsum += jetTrack->Pt();
 479       dR[iT] = recJets[iJ].DeltaR(jetTrack);
 480     }
 481
 482     fpHistPtAvE->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
 483     if(flagQ[iJ] == 1) fpHistPtAvEQ->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
 484     if(flagG[iJ] == 1) fpHistPtAvEG->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
 485
 486     if(iJ > 1) //fill mulit-jet histo
 487       fpHistPtAvE3->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
 488
 489     Int_t idxAOD[1000];
 490     TMath::Sort(nAODtracks, dR, idxAOD, kFALSE);
 491
 492     Double_t pTsum90Inv=0.;
 493     for(Int_t iT = 0; iT < nAODtracks; iT++){
 494       AliAODTrack * track = dynamic_cast<AliAODTrack*>(tra->At(idxAOD[iT]));
 495       if(!track) continue;
 496       pTsum90Inv += track->Pt();
 497       if(pTsum90Inv >= 0.9*pTsum){
 498         Double_t deltaR = recJets[iJ].DeltaR(track);
 499
 500         fpHistDrE->Fill(recJets[iJ].E(), deltaR, fWeight);
 501         if(flagQ[iJ] == 1) fpHistDrEQ->Fill(recJets[iJ].E(), deltaR, fWeight);
 502         if(flagG[iJ] == 1) fpHistDrEG->Fill(recJets[iJ].E(), deltaR, fWeight);
 503
 504         if(iJ > 1)
 505           fpHistDrE3->Fill(recJets[iJ].E(), deltaR, fWeight);
 506
 507         break;
 508       }
 509     }
 510   }
 511 }