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                                                             ,fAllRectotETDep(0)
  41                                                             ,fElectronMatchtotETDep(0)
  42                                                             ,fNeutralRectotET(0)
  43                                                             ,fTotEMRectotET(0)
  44
  45                                                             ,fMuonMatchtotETDep(0), fPionMatchtotETDep(0), fKaonMatchtotETDep(0), fProtonMatchtotETDep(0)
  46                                                             ,fTotChargedMatchtotETDep(0)
  47
  48                                                             ,fTotalRectotETDep(0)
  49
  50                                                             ,fESD(0)
  51 //                                                          ,fGeoUt(0)
  52
  53                                                             ,fHistAllRecETDep(0)
  54                                                             ,fHistAllRec(0)
  55                                                             ,fHistAllRectotETDep(0)
  56
  57                                                             ,fHistElectronRecETDep(0)
  58                                                             ,fHistElectronRec(0)
  59                                                             ,fHistElectronMatchtotETDep(0)
  60                                                             ,fHistElectronRecdEdxP(0)
  61
  62                                                             ,fHistNeutralRectotET(0)
  63
  64                                                             ,fHistTotEMRectotET(0)
  65
  66                                                             ,fHistMuonRecETDep(0)
  67                                                             ,fHistMuonRec(0)
  68                                                             ,fHistMuonMatchtotETDep(0)
  69                                                             ,fHistMuonRecdEdxP(0)
  70
  71                                                             ,fHistPionRecETDep(0)
  72                                                             ,fHistPionRec(0)
  73                                                             ,fHistPionMatchtotETDep(0)
  74                                                             ,fHistPionRecdEdxP(0)
  75
  76                                                             ,fHistKaonRecETDep(0)
  77                                                             ,fHistKaonRec(0)
  78                                                             ,fHistKaonMatchtotETDep(0)
  79                                                             ,fHistKaonRecdEdxP(0)
  80
  81                                                             ,fHistProtonRecETDep(0)
  82                                                             ,fHistProtonRec(0)
  83                                                             ,fHistProtonMatchtotETDep(0)
  84                                                             ,fHistProtonRecdEdxP(0)
  85
  86                                                             ,fHistTotChargedMatchtotETDep(0)
  87
  88                                                             ,fHistTotalRectotETDep(0)
  89
  90                                                             ,fHistDeltaRZ(0)
  91 {//constructor
  92   fHistogramNameSuffix = TString("EmcalRec");
  93
  94   fResCut = 0.02;
  95   //fResCut = fEmcalTrackDistanceCut;
  96
  97   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG));
  98   //TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  99   TGeoManager::Import("geometry.root");
 100 }
 101
 102 // dtor
 103 AliAnalysisEmEtReconstructed::~AliAnalysisEmEtReconstructed()
 104 {//Destructor
 105 //   delete fGeoUt;
 106
 107   delete fHistAllRecETDep;
 108   delete fHistAllRec;
 109   delete fHistAllRectotETDep;
 110
 111   delete fHistElectronRecETDep;
 112   delete fHistElectronRec;
 113   delete fHistElectronMatchtotETDep;
 114
 115   delete fHistElectronRecdEdxP;
 116
 117   delete fHistNeutralRectotET;
 118
 119   delete fHistTotEMRectotET;
 120
 121   delete fHistMuonRecETDep;
 122   delete fHistMuonRec;
 123   delete fHistMuonMatchtotETDep;
 124
 125   delete fHistMuonRecdEdxP;
 126
 127   delete fHistPionRecETDep;
 128   delete fHistPionRec;
 129   delete fHistPionMatchtotETDep;
 130
 131   delete fHistPionRecdEdxP;
 132
 133   delete fHistKaonRecETDep;
 134   delete fHistKaonRec;
 135   delete fHistKaonMatchtotETDep;
 136
 137   delete fHistKaonRecdEdxP;
 138
 139   delete fHistProtonRecETDep;
 140   delete fHistProtonRec;
 141   delete fHistProtonMatchtotETDep;
 142
 143   delete fHistProtonRecdEdxP;
 144
 145   delete fHistTotChargedMatchtotETDep;
 146
 147   delete fHistTotalRectotETDep;
 148
 149   //few checks
 150   delete fHistDeltaRZ;
 151
 152 }
 153
 154 Int_t AliAnalysisEmEtReconstructed::AnalyseEvent(AliVEvent* ev)
 155 { // analyse MC and real event info
 156   if(!ev){
 157     AliError("ERROR: Event does not exist");
 158     return 0;
 159   }
 160
 161   fESD = dynamic_cast<AliESDEvent*>(ev);
 162
 163 //   if(!fGeoUt){
 164 //     fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");//new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
 165 //     //fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 166 //     AliInfo("Creating new AliEMCALGeometry");
 167 //   }
 168 //   //fGeoUt = new AliEMCALGeometry("EMCAL_COMPLETE1","EMCAL");
 169 //   if(!fGeoUt){
 170 //     AliInfo("No fGeoUt!");
 171 //   }
 172 //   else{
 173 //     if(!fESD->GetEMCALMatrix(0)){
 174 //       AliInfo("No matrix!");
 175 //     }
 176 //     else{
 177 //       fGeoUt->SetMisalMatrix(fESD->GetEMCALMatrix(0),0);
 178 //     }
 179 //   }
 180
 181   fCentBin= -1;
 182   if(fDataSet==20100){//If this is Pb+Pb
 183     AliCentrality *centrality = ev->GetCentrality();
 184     if(fNCentBins<21) fCentBin= centrality->GetCentralityClass10(fCentralityMethod);
 185     else{ fCentBin= centrality->GetCentralityClass5(fCentralityMethod);}
 186   }
 187
 188   ResetEventValues();
 189
 190   // get all emcal clusters
 191   TRefArray* caloClusters = new TRefArray();
 192   fESD->GetEMCALClusters( caloClusters );
 193
 194   Int_t nCluster = caloClusters->GetEntries();
 195
 196   Float_t pos[3] = {0};
 197   TVector3 caloPos(0,0,0);
 198   TVector3 trackPos(0,0,0);
 199   Double_t res=0, delta_eta=0, delta_phi=0, maxPid=-99;
 200   Double_t xCluster[4]={0}, xCharged[7]={0};
 201
 202   AliESDtrack *track = 0;
 203
 204   // loop the clusters
 205   for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 206     {
 207       // Retrieve calo cluster information
 208       AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 209       Float_t caloE = caloCluster->E();
 210       caloCluster->GetPosition(pos);
 211       caloPos.SetXYZ(pos[0],pos[1],pos[2]);
 212
 213       // look for track that matches calo cluster
 214       //track = FindMatch(caloCluster, res); // Marcelo's matching
 215
 216       // *********************
 217       // tender's matching
 218       delta_eta = caloCluster->GetTrackDz();
 219       delta_phi = caloCluster->GetTrackDx();
 220
 221       if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
 222           track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
 223
 224       //if (track)
 225       //   if ( !fEsdtrackCutsITSTPC->IsSelected(track) )
 226       //       track = 0;
 227       // *********************
 228
 229       // Retrieve track PID
 230       if (track)
 231         maxPid = GetTrackPID(track);
 232       else
 233         maxPid = -99;
 234
 235       // calculate ET
 236       Double_t etDep = CalculateTransverseEnergy(caloCluster);
 237
 238       // All clusters
 239       //fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
 240       //fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
 241
 242       xCluster[0] = caloE;
 243       xCluster[1] = caloPos.Eta();
 244       xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
 245       xCluster[3] = caloCluster->GetNCells();
 246       fAllRectotETDep += etDep;
 247
 248
 249       if(fMakeSparse){
 250         fHistAllRecETDep->Fill(xCluster,etDep);
 251         fHistAllRec->Fill(xCluster);
 252       }
 253
 254       if (track)
 255         {
 256           xCharged[0] = track->Eta();
 257           xCharged[1] = track->Pt();
 258         }
 259       else
 260         {
 261           xCharged[0] = -99;
 262           xCharged[1] = -99;
 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           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           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           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           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           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 }