PWGLF/totEt/AliAnalysisEmEtReconstructed.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Base class for MC analysis
   4 //  - MC output
   5 //  implementation file
   6 //
   7 //*-- Author: Marcelo G. Munhoz (USP)
   8               //_________________________________________________________________________
   9
  10 #include "AliAnalysisEmEtReconstructed.h"
  11 #include "AliAnalysisEtCuts.h"
  12 #include "AliESDtrack.h"
  13 #include "AliStack.h"
  14 #include "AliVEvent.h"
  15 #include "AliMCEvent.h"
  16 #include "AliESDEvent.h"
  17 #include "TH2F.h"
  18 #include "TParticle.h"
  19 #include "AliGenHijingEventHeader.h"
  20 #include "AliGenPythiaEventHeader.h"
  21 #include "TList.h"
  22 #include "AliESDCaloCluster.h"
  23 #include "TGeoGlobalMagField.h"
  24 #include "AliMagF.h"
  25 #include "AliEMCALTrack.h"
  26 #include "AliESDtrackCuts.h"
  27 #include "AliEMCALGeometry.h"
  28 #include "AliExternalTrackParam.h"
  29 #include "AliTrackerBase.h"
  30 #include "TGeoManager.h"
  31 #include "AliCentrality.h"
  32
  33   using namespace std;
  34
  35 ClassImp(AliAnalysisEmEtReconstructed);
  36
  37
  38 // ctor
  39 AliAnalysisEmEtReconstructed::AliAnalysisEmEtReconstructed():AliAnalysisEtReconstructed()
  40                                                             ,fResCut(0)
  41                                                             ,fAllRectotETDep(0)
  42                                                             ,fElectronMatchtotETDep(0)
  43                                                             ,fNeutralRectotET(0)
  44                                                             ,fTotEMRectotET(0)
  45
  46                                                             ,fMuonMatchtotETDep(0), fPionMatchtotETDep(0), fKaonMatchtotETDep(0), fProtonMatchtotETDep(0)
  47                                                             ,fTotChargedMatchtotETDep(0)
  48
  49                                                             ,fTotalRectotETDep(0)
  50
  51                                                             ,fESD(0)
  52 //                                                          ,fGeoUt(0)
  53
  54                                                             ,fHistAllRecETDep(0)
  55                                                             ,fHistAllRec(0)
  56                                                             ,fHistAllRectotETDep(0)
  57
  58                                                             ,fHistElectronRecETDep(0)
  59                                                             ,fHistElectronRec(0)
  60                                                             ,fHistElectronMatchtotETDep(0)
  61                                                             ,fHistElectronRecdEdxP(0)
  62
  63                                                             ,fHistNeutralRectotET(0)
  64
  65                                                             ,fHistTotEMRectotET(0)
  66
  67                                                             ,fHistMuonRecETDep(0)
  68                                                             ,fHistMuonRec(0)
  69                                                             ,fHistMuonMatchtotETDep(0)
  70                                                             ,fHistMuonRecdEdxP(0)
  71
  72                                                             ,fHistPionRecETDep(0)
  73                                                             ,fHistPionRec(0)
  74                                                             ,fHistPionMatchtotETDep(0)
  75                                                             ,fHistPionRecdEdxP(0)
  76
  77                                                             ,fHistKaonRecETDep(0)
  78                                                             ,fHistKaonRec(0)
  79                                                             ,fHistKaonMatchtotETDep(0)
  80                                                             ,fHistKaonRecdEdxP(0)
  81
  82                                                             ,fHistProtonRecETDep(0)
  83                                                             ,fHistProtonRec(0)
  84                                                             ,fHistProtonMatchtotETDep(0)
  85                                                             ,fHistProtonRecdEdxP(0)
  86
  87                                                             ,fHistTotChargedMatchtotETDep(0)
  88
  89                                                             ,fHistTotalRectotETDep(0)
  90
  91                                                             ,fHistDeltaRZ(0)
  92 {//constructor
  93   fHistogramNameSuffix = TString("EmcalRec");
  94
  95   fResCut = 0.02;
  96   //fResCut = fEmcalTrackDistanceCut;
  97
  98   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG));
  99   //TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 100   TGeoManager::Import("geometry.root");
 101 }
 102
 103 // dtor
 104 AliAnalysisEmEtReconstructed::~AliAnalysisEmEtReconstructed()
 105 {//Destructor
 106 //   delete fGeoUt;
 107
 108   delete fHistAllRecETDep;
 109   delete fHistAllRec;
 110   delete fHistAllRectotETDep;
 111
 112   delete fHistElectronRecETDep;
 113   delete fHistElectronRec;
 114   delete fHistElectronMatchtotETDep;
 115
 116   delete fHistElectronRecdEdxP;
 117
 118   delete fHistNeutralRectotET;
 119
 120   delete fHistTotEMRectotET;
 121
 122   delete fHistMuonRecETDep;
 123   delete fHistMuonRec;
 124   delete fHistMuonMatchtotETDep;
 125
 126   delete fHistMuonRecdEdxP;
 127
 128   delete fHistPionRecETDep;
 129   delete fHistPionRec;
 130   delete fHistPionMatchtotETDep;
 131
 132   delete fHistPionRecdEdxP;
 133
 134   delete fHistKaonRecETDep;
 135   delete fHistKaonRec;
 136   delete fHistKaonMatchtotETDep;
 137
 138   delete fHistKaonRecdEdxP;
 139
 140   delete fHistProtonRecETDep;
 141   delete fHistProtonRec;
 142   delete fHistProtonMatchtotETDep;
 143
 144   delete fHistProtonRecdEdxP;
 145
 146   delete fHistTotChargedMatchtotETDep;
 147
 148   delete fHistTotalRectotETDep;
 149
 150   //few checks
 151   delete fHistDeltaRZ;
 152
 153 }
 154
 155 Int_t AliAnalysisEmEtReconstructed::AnalyseEvent(AliVEvent* ev)
 156 { // analyse MC and real event info
 157   if(!ev){
 158     AliError("ERROR: Event does not exist");
 159     return 0;
 160   }
 161
 162   fESD = dynamic_cast<AliESDEvent*>(ev);
 163   if (!fESD) {
 164     AliError("ERROR: Could not retrieve event");
 165     return 0;
 166   }
 167
 168 //   if(!fGeoUt){
 169 //     fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");//new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
 170 //     //fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 171 //     AliInfo("Creating new AliEMCALGeometry");
 172 //   }
 173 //   //fGeoUt = new AliEMCALGeometry("EMCAL_COMPLETE1","EMCAL");
 174 //   if(!fGeoUt){
 175 //     AliInfo("No fGeoUt!");
 176 //   }
 177 //   else{
 178 //     if(!fESD->GetEMCALMatrix(0)){
 179 //       AliInfo("No matrix!");
 180 //     }
 181 //     else{
 182 //       fGeoUt->SetMisalMatrix(fESD->GetEMCALMatrix(0),0);
 183 //     }
 184 //   }
 185
 186   fCentBin= -1;
 187   if(fDataSet==20100){//If this is Pb+Pb
 188     AliCentrality *centrality = ev->GetCentrality();
 189     if(fNCentBins<21) fCentBin= centrality->GetCentralityClass10(fCentralityMethod);
 190     else{ fCentBin= centrality->GetCentralityClass5(fCentralityMethod);}
 191   }
 192
 193   ResetEventValues();
 194
 195   // get all emcal clusters
 196   TRefArray* caloClusters = new TRefArray();
 197   fESD->GetEMCALClusters( caloClusters );
 198
 199   Int_t nCluster = caloClusters->GetEntries();
 200
 201   Float_t pos[3] = {0};
 202   TVector3 caloPos(0,0,0);
 203   TVector3 trackPos(0,0,0);
 204   Double_t res=0, delta_eta=0, delta_phi=0, maxPid=-99;
 205   Double_t xCluster[4]={0}, xCharged[7]={0};
 206
 207   AliESDtrack *track = 0;
 208
 209   // loop the clusters
 210   for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 211     {
 212       // Retrieve calo cluster information
 213       AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 214       Float_t caloE = caloCluster->E();
 215       caloCluster->GetPosition(pos);
 216       caloPos.SetXYZ(pos[0],pos[1],pos[2]);
 217
 218       // look for track that matches calo cluster
 219       //track = FindMatch(caloCluster, res); // Marcelo's matching
 220
 221       // *********************
 222       // tender's matching
 223       delta_eta = caloCluster->GetTrackDz();
 224       delta_phi = caloCluster->GetTrackDx();
 225
 226       if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
 227           track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
 228
 229       //if (track)
 230       //   if ( !fEsdtrackCutsITSTPC->IsSelected(track) )
 231       //       track = 0;
 232       // *********************
 233
 234       // Retrieve track PID
 235       if (track) {
 236         maxPid = GetTrackPID(track);
 237         xCharged[0] = track->Eta();
 238         xCharged[1] = track->Pt();      }
 239       else {
 240         maxPid = -99;
 241         xCharged[0] = -99;
 242         xCharged[1] = -99;
 243       }
 244
 245       // calculate ET
 246       Double_t etDep = CalculateTransverseEnergy(caloCluster);
 247
 248       // All clusters
 249       //fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
 250       //fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
 251
 252       xCluster[0] = caloE;
 253       xCluster[1] = caloPos.Eta();
 254       xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
 255       xCluster[3] = caloCluster->GetNCells();
 256       fAllRectotETDep += etDep;
 257
 258
 259       if(fMakeSparse){
 260         fHistAllRecETDep->Fill(xCluster,etDep);
 261         fHistAllRec->Fill(xCluster);
 262       }
 263
 264       xCharged[2] = caloE;
 265       xCharged[3] = caloPos.Eta();
 266       xCharged[4] = TMath::RadToDeg()*caloPos.Phi();
 267       xCharged[5] = caloCluster->GetNCells();
 268       xCharged[6] = res;
 269
 270       Bool_t isCharged = kFALSE;
 271
 272       if (maxPid == AliPID::kProton)
 273         {
 274
 275           if(fMakeSparse){
 276             fHistProtonRecETDep->Fill(xCharged,etDep);
 277             fHistProtonRec->Fill(xCharged);
 278           }
 279
 280           if (track) fHistProtonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 281
 282           if ((res>0.) && (res<fResCut))
 283             {
 284               fProtonMatchtotETDep += etDep;
 285
 286               isCharged = kTRUE;
 287             }
 288         }
 289       else if (maxPid == AliPID::kPion)
 290         {
 291
 292           if(fMakeSparse){
 293             fHistPionRecETDep->Fill(xCharged,etDep);
 294             fHistPionRec->Fill(xCharged);
 295           }
 296
 297           if (track) fHistPionRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 298
 299           if ((res>0.) && (res<fResCut))
 300             {
 301               fPionMatchtotETDep += etDep;
 302               isCharged = kTRUE;
 303             }
 304         }
 305       else if (maxPid == AliPID::kKaon)
 306         {
 307
 308           if(fMakeSparse){
 309             fHistKaonRecETDep->Fill(xCharged,etDep);
 310             fHistKaonRec->Fill(xCharged);
 311           }
 312
 313           if (track) fHistKaonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 314
 315           if ((res>0.) && (res<fResCut))
 316             {
 317
 318               fKaonMatchtotETDep += etDep;
 319               isCharged = kTRUE;
 320             }
 321         }
 322       else if (maxPid == AliPID::kMuon)
 323         {
 324
 325           if(fMakeSparse){
 326             fHistMuonRecETDep->Fill(xCharged,etDep);
 327             fHistMuonRec->Fill(xCharged);
 328           }
 329           if (track) fHistMuonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 330
 331           if ((res>0.) && (res<fResCut))
 332             {
 333
 334               fMuonMatchtotETDep += etDep;
 335               isCharged = kTRUE;
 336             }
 337         }
 338       else if (maxPid == AliPID::kElectron)
 339         {
 340
 341           if(fMakeSparse){
 342             fHistElectronRecETDep->Fill(xCharged,etDep);
 343             fHistElectronRec->Fill(xCharged);
 344           }
 345
 346           if (track) fHistElectronRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 347
 348           if ((res>0.) && (res<fResCut))
 349             {
 350               fElectronMatchtotETDep += etDep;
 351               isCharged = kTRUE;
 352             }
 353         }
 354
 355       if (!isCharged)
 356         {
 357           fNeutralRectotET += etDep;
 358         }
 359
 360     } // end of loop over clusters
 361
 362   fTotEMRectotET = fElectronMatchtotETDep + fNeutralRectotET;
 363   fTotChargedMatchtotETDep = fMuonMatchtotETDep + fPionMatchtotETDep + fKaonMatchtotETDep + fProtonMatchtotETDep;
 364   fTotalRectotETDep = fTotEMRectotET + fTotChargedMatchtotETDep;
 365
 366   fHistAllRectotETDep->Fill(fAllRectotETDep);
 367
 368   fHistElectronMatchtotETDep->Fill(fElectronMatchtotETDep);
 369   fHistNeutralRectotET->Fill(fNeutralRectotET);
 370
 371   fHistTotEMRectotET->Fill(fTotEMRectotET);
 372
 373   fHistMuonMatchtotETDep->Fill(fMuonMatchtotETDep);
 374   fHistPionMatchtotETDep->Fill(fPionMatchtotETDep);
 375   fHistKaonMatchtotETDep->Fill(fKaonMatchtotETDep);
 376   fHistProtonMatchtotETDep->Fill(fProtonMatchtotETDep);
 377   fHistTotChargedMatchtotETDep->Fill(fTotChargedMatchtotETDep);
 378
 379   fHistTotalRectotETDep->Fill(fTotalRectotETDep);
 380
 381   delete caloClusters;
 382
 383   return 0;
 384 }
 385
 386 void AliAnalysisEmEtReconstructed::Init()
 387 { // init
 388   AliAnalysisEt::Init();
 389 }
 390
 391
 392 void AliAnalysisEmEtReconstructed::ResetEventValues()
 393 { // reset event values
 394   AliAnalysisEt::ResetEventValues();
 395
 396   // collision geometry defaults for p+p:
 397   fAllRectotETDep = 0;
 398
 399   fElectronMatchtotETDep = 0;
 400   fNeutralRectotET = 0;
 401
 402   fTotEMRectotET = 0;
 403
 404   fMuonMatchtotETDep = 0; fPionMatchtotETDep = 0; fKaonMatchtotETDep = 0; fProtonMatchtotETDep = 0;
 405   fTotChargedMatchtotETDep = 0;
 406
 407   fTotalRectotETDep = 0;
 408 }
 409
 410
 411 void AliAnalysisEmEtReconstructed::CreateHistograms()
 412 { // histogram related additions
 413   //AliAnalysisEt::CreateHistograms();
 414
 415   if(fMakeSparse){
 416     fHistAllRecETDep = CreateClusterHistoSparse("fHistAllRecETDep_","E_{T}, all particles");
 417     fHistAllRec = CreateClusterHistoSparse("fHistAllRec_","counts, all particles");
 418   }
 419   TString histname = "fHistAllRectotETDep_" + fHistogramNameSuffix;
 420   fHistAllRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
 421
 422
 423   if(fMakeSparse){
 424     fHistElectronRecETDep = CreateChargedPartHistoSparse("fHistElectronRecETDep_","E_{T}, electrons");
 425     fHistElectronRec = CreateChargedPartHistoSparse("fHistElectronRec_","counts, electrons");
 426   }
 427   histname = "fHistElectronMatchtotETDep_" + fHistogramNameSuffix;
 428   fHistElectronMatchtotETDep = new TH1F(histname.Data(),"total ET, MC primary Electrons",fgNumOfEBins, fgEAxis);
 429
 430   histname = "fHistElectronRecdEdxP_" + fHistogramNameSuffix;
 431   fHistElectronRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 432
 433   histname = "fHistNeutralRectotET_" + fHistogramNameSuffix;
 434   fHistNeutralRectotET = new TH1F(histname.Data(),"total ET, neutral particles",fgNumOfEBins, fgEAxis);
 435
 436   histname = "fHistTotEMRectotET_" + fHistogramNameSuffix;
 437   fHistTotEMRectotET = new TH1F(histname.Data(),"total electromagnetic ET",fgNumOfEBins, fgEAxis);
 438
 439   if(fMakeSparse){
 440     fHistMuonRecETDep = CreateChargedPartHistoSparse("fHistMuonRecETDep_","E_{T}, muons");
 441     fHistMuonRec = CreateChargedPartHistoSparse("fHistMuonRec_","counts, muons");
 442   }
 443   histname = "fHistMuonMatchtotETDep_" + fHistogramNameSuffix;
 444   fHistMuonMatchtotETDep = new TH1F(histname.Data(),"total ET, Muons",fgNumOfEBins, fgEAxis);
 445
 446   histname = "fHistMuonRecdEdxP_" + fHistogramNameSuffix;
 447   fHistMuonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 448
 449   if(fMakeSparse){
 450     fHistPionRecETDep = CreateChargedPartHistoSparse("fHistPionRecETDep_","E_{T}, pions");
 451     fHistPionRec = CreateChargedPartHistoSparse("fHistPionRec_","counts, pions");
 452   }
 453   histname = "fHistPionMatchtotETDep_" + fHistogramNameSuffix;
 454   fHistPionMatchtotETDep = new TH1F(histname.Data(),"total ET, Pions",fgNumOfEBins, fgEAxis);
 455   histname = "fHistPionRecdEdxP_" + fHistogramNameSuffix;
 456   fHistPionRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 457
 458   if(fMakeSparse){
 459     fHistKaonRecETDep = CreateChargedPartHistoSparse("fHistKaonRecETDep_","E_{T}, kaons");
 460     fHistKaonRec = CreateChargedPartHistoSparse("fHistKaonRec_","counts, kaons");
 461   }
 462   histname = "fHistKaonMatchtotETDep_" + fHistogramNameSuffix;
 463   fHistKaonMatchtotETDep = new TH1F(histname.Data(),"total ET, Kaons",fgNumOfEBins, fgEAxis);
 464
 465   histname = "fHistKaonRecdEdxP_" + fHistogramNameSuffix;
 466   fHistKaonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 467
 468   if(fMakeSparse){
 469     fHistProtonRecETDep = CreateChargedPartHistoSparse("fHistProtonRecETDep_","E_{T}, protons");
 470     fHistProtonRec = CreateChargedPartHistoSparse("fHistProtonRec_","counts, protons");
 471   }
 472   histname = "fHistProtonMatchtotETDep_" + fHistogramNameSuffix;
 473   fHistProtonMatchtotETDep = new TH1F(histname.Data(),"total ET, Protons",fgNumOfEBins, fgEAxis);
 474
 475   histname = "fHistProtonRecdEdxP_" + fHistogramNameSuffix;
 476   fHistProtonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 477
 478   histname = "fHistTotChargedMatchtotETDep_" + fHistogramNameSuffix;
 479   fHistTotChargedMatchtotETDep = new TH1F(histname.Data(),"total ET, charged particles",fgNumOfEBins, fgEAxis);
 480
 481   histname = "fHistTotalRectotETDep_" + fHistogramNameSuffix;
 482   fHistTotalRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
 483
 484   histname = "fHistDeltaRZ_" + fHistogramNameSuffix;
 485   fHistDeltaRZ = new TH2F(histname,"#Delta#phi vs #Delta#eta (track projection - cluster position)",200,-0.1,0.1,200,-0.1,0.1);
 486 }
 487
 488 void AliAnalysisEmEtReconstructed::FillOutputList(TList *list)
 489 {//Function for filling the output list
 490   //AliAnalysisEt::FillOutputList(list);
 491
 492
 493   if(fMakeSparse){
 494     list->Add(fHistAllRecETDep);
 495     list->Add(fHistAllRec);
 496   }
 497   list->Add(fHistAllRectotETDep);
 498
 499   if(fMakeSparse){
 500     list->Add(fHistElectronRecETDep);
 501     list->Add(fHistElectronRec);
 502   }
 503   list->Add(fHistElectronMatchtotETDep);
 504   list->Add(fHistElectronRecdEdxP);
 505
 506
 507   list->Add(fHistTotEMRectotET);
 508
 509   list->Add(fHistMuonRec);
 510   list->Add(fHistMuonRecdEdxP);
 511   list->Add(fHistMuonMatchtotETDep);
 512
 513   if(fMakeSparse){
 514     list->Add(fHistPionRecETDep);
 515     list->Add(fHistPionRec);
 516   }
 517   list->Add(fHistPionMatchtotETDep);
 518   list->Add(fHistPionRecdEdxP);
 519
 520   if(fMakeSparse){
 521     list->Add(fHistKaonRecETDep);
 522     list->Add(fHistKaonRec);
 523   }
 524   list->Add(fHistKaonMatchtotETDep);
 525   list->Add(fHistKaonRecdEdxP);
 526
 527   if(fMakeSparse){
 528     list->Add(fHistProtonRecETDep);
 529     list->Add(fHistProtonRec);
 530   }
 531   list->Add(fHistProtonMatchtotETDep);
 532   list->Add(fHistProtonRecdEdxP);
 533
 534   list->Add(fHistTotChargedMatchtotETDep);
 535   list->Add(fHistTotalRectotETDep);
 536
 537   list->Add(fHistDeltaRZ);
 538 }
 539
 540 //________________________________________________________________________
 541 Double_t AliAnalysisEmEtReconstructed::GetTrackPID(const AliESDtrack *track) const
 542 {//Get the default track ID
 543   const Double_t *pidWeights = track->PID();
 544   Int_t maxpid = -1;
 545   Double_t maxpidweight = 0;
 546
 547   if (pidWeights)
 548     {
 549       for (Int_t p =0; p < AliPID::kSPECIES; p++)
 550         {
 551           if (pidWeights[p] > maxpidweight)
 552             {
 553               maxpidweight = pidWeights[p];
 554               maxpid = p;
 555             }
 556         }
 557     }
 558
 559   return maxpid;
 560 }