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 res=0, maxPid=-99;
 206   Double_t xCluster[4]={0}, xCharged[7]={0};
 207
 208   AliESDtrack *track = 0;
 209
 210   // loop the clusters
 211   for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 212     {
 213       // Retrieve calo cluster information
 214       AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 215       Float_t caloE = caloCluster->E();
 216       caloCluster->GetPosition(pos);
 217       caloPos.SetXYZ(pos[0],pos[1],pos[2]);
 218
 219       // look for track that matches calo cluster
 220       //track = FindMatch(caloCluster, res); // Marcelo's matching
 221
 222       // *********************
 223       // tender's matching
 224 //       delta_eta = caloCluster->GetTrackDz();
 225 //       delta_phi = caloCluster->GetTrackDx();
 226
 227       if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
 228           track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
 229
 230       //if (track)
 231       //   if ( !fEsdtrackCutsITSTPC->IsSelected(track) )
 232       //       track = 0;
 233       // *********************
 234
 235       // Retrieve track PID
 236       if (track) {
 237         maxPid = GetTrackPID(track);
 238         xCharged[0] = track->Eta();
 239         xCharged[1] = track->Pt();      }
 240       else {
 241         maxPid = -99;
 242         xCharged[0] = -99;
 243         xCharged[1] = -99;
 244       }
 245
 246       // calculate ET
 247       Double_t etDep = CorrectForReconstructionEfficiency(*caloCluster);
 248
 249       // All clusters
 250       //fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
 251       //fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
 252
 253       xCluster[0] = caloE;
 254       xCluster[1] = caloPos.Eta();
 255       xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
 256       xCluster[3] = caloCluster->GetNCells();
 257       fAllRectotETDep += etDep;
 258
 259
 260       if(fMakeSparse){
 261         fHistAllRecETDep->Fill(xCluster,etDep);
 262         fHistAllRec->Fill(xCluster);
 263       }
 264
 265       xCharged[2] = caloE;
 266       xCharged[3] = caloPos.Eta();
 267       xCharged[4] = TMath::RadToDeg()*caloPos.Phi();
 268       xCharged[5] = caloCluster->GetNCells();
 269       xCharged[6] = res;
 270
 271       Bool_t isCharged = kFALSE;
 272
 273       if (maxPid == AliPID::kProton)
 274         {
 275
 276           if(fMakeSparse){
 277             fHistProtonRecETDep->Fill(xCharged,etDep);
 278             fHistProtonRec->Fill(xCharged);
 279           }
 280
 281           if (track) fHistProtonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 282
 283           if ((res>0.) && (res<fResCut))
 284             {
 285               fProtonMatchtotETDep += etDep;
 286
 287               isCharged = kTRUE;
 288             }
 289         }
 290       else if (maxPid == AliPID::kPion)
 291         {
 292
 293           if(fMakeSparse){
 294             fHistPionRecETDep->Fill(xCharged,etDep);
 295             fHistPionRec->Fill(xCharged);
 296           }
 297
 298           if (track) fHistPionRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 299
 300           if ((res>0.) && (res<fResCut))
 301             {
 302               fPionMatchtotETDep += etDep;
 303               isCharged = kTRUE;
 304             }
 305         }
 306       else if (maxPid == AliPID::kKaon)
 307         {
 308
 309           if(fMakeSparse){
 310             fHistKaonRecETDep->Fill(xCharged,etDep);
 311             fHistKaonRec->Fill(xCharged);
 312           }
 313
 314           if (track) fHistKaonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 315
 316           if ((res>0.) && (res<fResCut))
 317             {
 318
 319               fKaonMatchtotETDep += etDep;
 320               isCharged = kTRUE;
 321             }
 322         }
 323       else if (maxPid == AliPID::kMuon)
 324         {
 325
 326           if(fMakeSparse){
 327             fHistMuonRecETDep->Fill(xCharged,etDep);
 328             fHistMuonRec->Fill(xCharged);
 329           }
 330           if (track) fHistMuonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 331
 332           if ((res>0.) && (res<fResCut))
 333             {
 334
 335               fMuonMatchtotETDep += etDep;
 336               isCharged = kTRUE;
 337             }
 338         }
 339       else if (maxPid == AliPID::kElectron)
 340         {
 341
 342           if(fMakeSparse){
 343             fHistElectronRecETDep->Fill(xCharged,etDep);
 344             fHistElectronRec->Fill(xCharged);
 345           }
 346
 347           if (track) fHistElectronRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 348
 349           if ((res>0.) && (res<fResCut))
 350             {
 351               fElectronMatchtotETDep += etDep;
 352               isCharged = kTRUE;
 353             }
 354         }
 355
 356       if (!isCharged)
 357         {
 358           fNeutralRectotET += etDep;
 359         }
 360
 361     } // end of loop over clusters
 362
 363   fTotEMRectotET = fElectronMatchtotETDep + fNeutralRectotET;
 364   fTotChargedMatchtotETDep = fMuonMatchtotETDep + fPionMatchtotETDep + fKaonMatchtotETDep + fProtonMatchtotETDep;
 365   fTotalRectotETDep = fTotEMRectotET + fTotChargedMatchtotETDep;
 366
 367   fHistAllRectotETDep->Fill(fAllRectotETDep);
 368
 369   fHistElectronMatchtotETDep->Fill(fElectronMatchtotETDep);
 370   fHistNeutralRectotET->Fill(fNeutralRectotET);
 371
 372   fHistTotEMRectotET->Fill(fTotEMRectotET);
 373
 374   fHistMuonMatchtotETDep->Fill(fMuonMatchtotETDep);
 375   fHistPionMatchtotETDep->Fill(fPionMatchtotETDep);
 376   fHistKaonMatchtotETDep->Fill(fKaonMatchtotETDep);
 377   fHistProtonMatchtotETDep->Fill(fProtonMatchtotETDep);
 378   fHistTotChargedMatchtotETDep->Fill(fTotChargedMatchtotETDep);
 379
 380   fHistTotalRectotETDep->Fill(fTotalRectotETDep);
 381
 382   delete caloClusters;
 383
 384   return 0;
 385 }
 386
 387 void AliAnalysisEmEtReconstructed::Init()
 388 { // init
 389   AliAnalysisEt::Init();
 390 }
 391
 392
 393 void AliAnalysisEmEtReconstructed::ResetEventValues()
 394 { // reset event values
 395   AliAnalysisEt::ResetEventValues();
 396
 397   // collision geometry defaults for p+p:
 398   fAllRectotETDep = 0;
 399
 400   fElectronMatchtotETDep = 0;
 401   fNeutralRectotET = 0;
 402
 403   fTotEMRectotET = 0;
 404
 405   fMuonMatchtotETDep = 0; fPionMatchtotETDep = 0; fKaonMatchtotETDep = 0; fProtonMatchtotETDep = 0;
 406   fTotChargedMatchtotETDep = 0;
 407
 408   fTotalRectotETDep = 0;
 409 }
 410
 411
 412 void AliAnalysisEmEtReconstructed::CreateHistograms()
 413 { // histogram related additions
 414   //AliAnalysisEt::CreateHistograms();
 415
 416   if(fMakeSparse){
 417     fHistAllRecETDep = CreateClusterHistoSparse("fHistAllRecETDep_","E_{T}, all particles");
 418     fHistAllRec = CreateClusterHistoSparse("fHistAllRec_","counts, all particles");
 419   }
 420   TString histname = "fHistAllRectotETDep_" + fHistogramNameSuffix;
 421   fHistAllRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
 422
 423
 424   if(fMakeSparse){
 425     fHistElectronRecETDep = CreateChargedPartHistoSparse("fHistElectronRecETDep_","E_{T}, electrons");
 426     fHistElectronRec = CreateChargedPartHistoSparse("fHistElectronRec_","counts, electrons");
 427   }
 428   histname = "fHistElectronMatchtotETDep_" + fHistogramNameSuffix;
 429   fHistElectronMatchtotETDep = new TH1F(histname.Data(),"total ET, MC primary Electrons",fgNumOfEBins, fgEAxis);
 430
 431   histname = "fHistElectronRecdEdxP_" + fHistogramNameSuffix;
 432   fHistElectronRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 433
 434   histname = "fHistNeutralRectotET_" + fHistogramNameSuffix;
 435   fHistNeutralRectotET = new TH1F(histname.Data(),"total ET, neutral particles",fgNumOfEBins, fgEAxis);
 436
 437   histname = "fHistTotEMRectotET_" + fHistogramNameSuffix;
 438   fHistTotEMRectotET = new TH1F(histname.Data(),"total electromagnetic ET",fgNumOfEBins, fgEAxis);
 439
 440   if(fMakeSparse){
 441     fHistMuonRecETDep = CreateChargedPartHistoSparse("fHistMuonRecETDep_","E_{T}, muons");
 442     fHistMuonRec = CreateChargedPartHistoSparse("fHistMuonRec_","counts, muons");
 443   }
 444   histname = "fHistMuonMatchtotETDep_" + fHistogramNameSuffix;
 445   fHistMuonMatchtotETDep = new TH1F(histname.Data(),"total ET, Muons",fgNumOfEBins, fgEAxis);
 446
 447   histname = "fHistMuonRecdEdxP_" + fHistogramNameSuffix;
 448   fHistMuonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 449
 450   if(fMakeSparse){
 451     fHistPionRecETDep = CreateChargedPartHistoSparse("fHistPionRecETDep_","E_{T}, pions");
 452     fHistPionRec = CreateChargedPartHistoSparse("fHistPionRec_","counts, pions");
 453   }
 454   histname = "fHistPionMatchtotETDep_" + fHistogramNameSuffix;
 455   fHistPionMatchtotETDep = new TH1F(histname.Data(),"total ET, Pions",fgNumOfEBins, fgEAxis);
 456   histname = "fHistPionRecdEdxP_" + fHistogramNameSuffix;
 457   fHistPionRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 458
 459   if(fMakeSparse){
 460     fHistKaonRecETDep = CreateChargedPartHistoSparse("fHistKaonRecETDep_","E_{T}, kaons");
 461     fHistKaonRec = CreateChargedPartHistoSparse("fHistKaonRec_","counts, kaons");
 462   }
 463   histname = "fHistKaonMatchtotETDep_" + fHistogramNameSuffix;
 464   fHistKaonMatchtotETDep = new TH1F(histname.Data(),"total ET, Kaons",fgNumOfEBins, fgEAxis);
 465
 466   histname = "fHistKaonRecdEdxP_" + fHistogramNameSuffix;
 467   fHistKaonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 468
 469   if(fMakeSparse){
 470     fHistProtonRecETDep = CreateChargedPartHistoSparse("fHistProtonRecETDep_","E_{T}, protons");
 471     fHistProtonRec = CreateChargedPartHistoSparse("fHistProtonRec_","counts, protons");
 472   }
 473   histname = "fHistProtonMatchtotETDep_" + fHistogramNameSuffix;
 474   fHistProtonMatchtotETDep = new TH1F(histname.Data(),"total ET, Protons",fgNumOfEBins, fgEAxis);
 475
 476   histname = "fHistProtonRecdEdxP_" + fHistogramNameSuffix;
 477   fHistProtonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 478
 479   histname = "fHistTotChargedMatchtotETDep_" + fHistogramNameSuffix;
 480   fHistTotChargedMatchtotETDep = new TH1F(histname.Data(),"total ET, charged particles",fgNumOfEBins, fgEAxis);
 481
 482   histname = "fHistTotalRectotETDep_" + fHistogramNameSuffix;
 483   fHistTotalRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
 484
 485   histname = "fHistDeltaRZ_" + fHistogramNameSuffix;
 486   fHistDeltaRZ = new TH2F(histname,"#Delta#phi vs #Delta#eta (track projection - cluster position)",200,-0.1,0.1,200,-0.1,0.1);
 487 }
 488
 489 void AliAnalysisEmEtReconstructed::FillOutputList(TList *list)
 490 {//Function for filling the output list
 491   //AliAnalysisEt::FillOutputList(list);
 492
 493
 494   if(fMakeSparse){
 495     list->Add(fHistAllRecETDep);
 496     list->Add(fHistAllRec);
 497   }
 498   list->Add(fHistAllRectotETDep);
 499
 500   if(fMakeSparse){
 501     list->Add(fHistElectronRecETDep);
 502     list->Add(fHistElectronRec);
 503   }
 504   list->Add(fHistElectronMatchtotETDep);
 505   list->Add(fHistElectronRecdEdxP);
 506
 507
 508   list->Add(fHistTotEMRectotET);
 509
 510   list->Add(fHistMuonRec);
 511   list->Add(fHistMuonRecdEdxP);
 512   list->Add(fHistMuonMatchtotETDep);
 513
 514   if(fMakeSparse){
 515     list->Add(fHistPionRecETDep);
 516     list->Add(fHistPionRec);
 517   }
 518   list->Add(fHistPionMatchtotETDep);
 519   list->Add(fHistPionRecdEdxP);
 520
 521   if(fMakeSparse){
 522     list->Add(fHistKaonRecETDep);
 523     list->Add(fHistKaonRec);
 524   }
 525   list->Add(fHistKaonMatchtotETDep);
 526   list->Add(fHistKaonRecdEdxP);
 527
 528   if(fMakeSparse){
 529     list->Add(fHistProtonRecETDep);
 530     list->Add(fHistProtonRec);
 531   }
 532   list->Add(fHistProtonMatchtotETDep);
 533   list->Add(fHistProtonRecdEdxP);
 534
 535   list->Add(fHistTotChargedMatchtotETDep);
 536   list->Add(fHistTotalRectotETDep);
 537
 538   list->Add(fHistDeltaRZ);
 539 }
 540
 541 //________________________________________________________________________
 542 Double_t AliAnalysisEmEtReconstructed::GetTrackPID(const AliESDtrack *track) const
 543 {//Get the default track ID
 544   const Double_t *pidWeights = track->PID();
 545   Int_t maxpid = -1;
 546   Double_t maxpidweight = 0;
 547
 548   if (pidWeights)
 549     {
 550       for (Int_t p =0; p < AliPID::kSPECIES; p++)
 551         {
 552           if (pidWeights[p] > maxpidweight)
 553             {
 554               maxpidweight = pidWeights[p];
 555               maxpid = p;
 556             }
 557         }
 558     }
 559
 560   return maxpid;
 561 }