PWG4/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
  32 using namespace std;
  33
  34 ClassImp(AliAnalysisEmEtReconstructed);
  35
  36
  37 // ctor
  38 AliAnalysisEmEtReconstructed::AliAnalysisEmEtReconstructed():AliAnalysisEtReconstructed()
  39 ,fAllRectotETDep(0)
  40 ,fElectronMatchtotETDep(0)
  41 ,fNeutralRectotET(0)
  42 ,fTotEMRectotET(0)
  43
  44 ,fMuonMatchtotETDep(0), fPionMatchtotETDep(0), fKaonMatchtotETDep(0), fProtonMatchtotETDep(0)
  45 ,fTotChargedMatchtotETDep(0)
  46
  47 ,fTotalRectotETDep(0)
  48
  49 ,fESD(0)
  50 ,fGeoUt(0)
  51
  52 //,fHistAllRecEtaEDepETDep(0)
  53 //,fHistAllRecEtaETDep(0)
  54 ,fHistAllRecETDep(0)
  55 ,fHistAllRec(0)
  56 ,fHistAllRectotETDep(0)
  57
  58 /*
  59 ,fHistElectronMatchEtaEDepETDep(0)
  60 ,fHistElectronMatchEtaPtETDep(0)
  61 ,fHistElectronMatchEtaETDep(0)
  62 ,fHistElectronMatchEtaPt(0)
  63
  64 ,fHistElectronRec_ResEDep_ETDep(0)
  65 ,fHistElectronRec_ResPt_ETDep(0)
  66 ,fHistElectronRec_ResEDep(0)
  67 ,fHistElectronRec_ResPt(0)
  68 */
  69 ,fHistElectronRecETDep(0)
  70 ,fHistElectronRec(0)
  71 ,fHistElectronMatchtotETDep(0)
  72 ,fHistElectronRecdEdxP(0)
  73
  74 /*
  75 ,fHistNeutralRec_EtaE_ET(0)
  76 //,fHistNeutralRec_EtaPt_ET(0)
  77 ,fHistNeutralRec_EtaET(0)
  78 ,fHistNeutralRec_EtaE(0)
  79 //,fHistNeutralRec_EtaPt(0)
  80 */
  81 ,fHistNeutralRectotET(0)
  82
  83 ,fHistTotEMRectotET(0)
  84
  85 /*
  86 ,fHistMuonMatchEtaEDepETDep(0)
  87 ,fHistMuonMatchEtaPtETDep(0)
  88 ,fHistMuonMatchEtaETDep(0)
  89 ,fHistMuonMatchEtaPt(0)
  90
  91 ,fHistMuonRecResEDepETDep(0)
  92 ,fHistMuonRecResPtETDep(0)
  93 ,fHistMuonRecResEDep(0)
  94 ,fHistMuonRecResPt(0)
  95 */
  96 ,fHistMuonRecETDep(0)
  97 ,fHistMuonRec(0)
  98 ,fHistMuonMatchtotETDep(0)
  99 ,fHistMuonRecdEdxP(0)
 100
 101 /*
 102 ,fHistPionMatchEtaEDepETDep(0)
 103 ,fHistPionMatchEtaPtETDep(0)
 104 ,fHistPionMatchEtaETDep(0)
 105 ,fHistPionMatchEtaPt(0)
 106
 107 ,fHistPionRecResEDepETDep(0)
 108 ,fHistPionRecResPtETDep(0)
 109 ,fHistPionRecResEDep(0)
 110 ,fHistPionRecResPt(0)
 111 */
 112 ,fHistPionRecETDep(0)
 113 ,fHistPionRec(0)
 114 ,fHistPionMatchtotETDep(0)
 115 ,fHistPionRecdEdxP(0)
 116
 117 /*
 118 ,fHistKaonMatchEtaEDepETDep(0)
 119 ,fHistKaonMatchEtaPtETDep(0)
 120 ,fHistKaonMatchEtaETDep(0)
 121 ,fHistKaonMatchEtaPt(0)
 122
 123 ,fHistKaonRecResEDepETDep(0)
 124 ,fHistKaonRecResPtETDep(0)
 125 ,fHistKaonRecResEDep(0)
 126 ,fHistKaonRecResPt(0)
 127 */
 128 ,fHistKaonRecETDep(0)
 129 ,fHistKaonRec(0)
 130 ,fHistKaonMatchtotETDep(0)
 131 ,fHistKaonRecdEdxP(0)
 132
 133 /*
 134 ,fHistProtonMatchEtaEDepETDep(0)
 135 ,fHistProtonMatchEtaPtETDep(0)
 136 ,fHistProtonMatchEtaETDep(0)
 137 ,fHistProtonMatchEtaPt(0)
 138
 139 ,fHistProtonRecResEDepETDep(0)
 140 ,fHistProtonRecResPtETDep(0)
 141 ,fHistProtonRecResEDep(0)
 142 ,fHistProtonRecResPt(0)
 143 */
 144 ,fHistProtonRecETDep(0)
 145 ,fHistProtonRec(0)
 146 ,fHistProtonMatchtotETDep(0)
 147 ,fHistProtonRecdEdxP(0)
 148
 149 ,fHistTotChargedMatchtotETDep(0)
 150
 151 ,fHistTotalRectotETDep(0)
 152
 153 ,fHistDeltaRZ(0)
 154 {//constructor
 155         fHistogramNameSuffix = TString("EmcalRec");
 156
 157         fResCut = 0.02;
 158         //fResCut = fEmcalTrackDistanceCut;
 159
 160         TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG));
 161         //TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 162         TGeoManager::Import("geometry.root");
 163         //fGeoUt = new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
 164 }
 165
 166 // dtor
 167 AliAnalysisEmEtReconstructed::~AliAnalysisEmEtReconstructed()
 168 {//Destructor
 169   //Marcelo, are you sure you clean up all memory?
 170   delete fGeoUt;
 171
 172         delete fHistAllRecETDep;//Marcelo please add comment
 173         delete fHistAllRec;//Marcelo please add comment
 174         delete fHistAllRectotETDep;//Marcelo please add comment
 175
 176         delete fHistElectronRecETDep;//Marcelo please add comment
 177         delete fHistElectronRec;//Marcelo please add comment
 178         delete fHistElectronMatchtotETDep;//Marcelo please add comment
 179
 180         delete fHistElectronRecdEdxP;//Marcelo please add comment
 181
 182         delete fHistNeutralRectotET;//Marcelo please add comment
 183
 184         delete fHistTotEMRectotET;//Marcelo please add comment
 185
 186         delete fHistMuonRecETDep;//Marcelo please add comment
 187         delete fHistMuonRec;//Marcelo please add comment
 188         delete fHistMuonMatchtotETDep;//Marcelo please add comment
 189
 190         delete fHistMuonRecdEdxP;//Marcelo please add comment
 191
 192         delete fHistPionRecETDep;//Marcelo please add comment
 193         delete fHistPionRec;//Marcelo please add comment
 194         delete fHistPionMatchtotETDep;//Marcelo please add comment
 195
 196         delete fHistPionRecdEdxP;//Marcelo please add comment
 197
 198         delete fHistKaonRecETDep;//Marcelo please add comment
 199         delete fHistKaonRec;//Marcelo please add comment
 200         delete fHistKaonMatchtotETDep;//Marcelo please add comment
 201
 202         delete fHistKaonRecdEdxP;//Marcelo please add comment
 203
 204         delete fHistProtonRecETDep;//Marcelo please add comment
 205         delete fHistProtonRec;//Marcelo please add comment
 206         delete fHistProtonMatchtotETDep;//Marcelo please add comment
 207
 208         delete fHistProtonRecdEdxP;//Marcelo please add comment
 209
 210         delete fHistTotChargedMatchtotETDep;//Marcelo please add comment
 211
 212         delete fHistTotalRectotETDep;//Marcelo please add comment
 213
 214         //few checks
 215         delete fHistDeltaRZ;//Marcelo please add comment
 216         //TH2F * TH1F * delete
 217
 218 }
 219
 220 Int_t AliAnalysisEmEtReconstructed::AnalyseEvent(AliVEvent* ev)
 221 { // analyse MC and real event info
 222         if(!ev){
 223                 Printf("ERROR: Event does not exist");
 224                 return 0;
 225         }
 226
 227         fESD = dynamic_cast<AliESDEvent*>(ev);
 228
 229         fGeoUt = new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
 230         //fGeoUt = new AliEMCALGeometry("EMCAL_COMPLETE1","EMCAL");
 231         fGeoUt->SetMisalMatrix(fESD->GetEMCALMatrix(0),0);
 232
 233         ResetEventValues();
 234
 235         // get all emcal clusters
 236         TRefArray* caloClusters = new TRefArray();
 237         fESD->GetEMCALClusters( caloClusters );
 238
 239         Int_t nCluster = caloClusters->GetEntries();
 240
 241         Float_t pos[3] = {0};
 242         TVector3 caloPos(0,0,0);
 243         TVector3 trackPos(0,0,0);
 244         Double_t res=0, maxPid=-99;
 245         Double_t xCluster[4]={0}, xCharged[7]={0};
 246
 247         // loop the clusters
 248         for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 249         {
 250                 // Retrieve calo cluster information
 251                 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 252                 Float_t caloE = caloCluster->E();
 253                 caloCluster->GetPosition(pos);
 254                 caloPos.SetXYZ(pos[0],pos[1],pos[2]);
 255
 256                 // look for track that matches calo cluster
 257                 AliESDtrack *track = FindMatch(caloCluster, res);
 258
 259                 // Retrieve track PID
 260                 if (track)
 261                         maxPid = GetTrackPID(track);
 262                 else
 263                         maxPid = -99;
 264
 265                 // calculate ET
 266                 Double_t etDep = CalculateTransverseEnergy(caloCluster);
 267
 268                 // All clusters
 269                 //fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
 270                 //fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
 271
 272                 xCluster[0] = caloE;
 273                 xCluster[1] = caloPos.Eta();
 274                 xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
 275                 xCluster[3] = caloCluster->GetNCells();
 276                 fAllRectotETDep += etDep;
 277
 278                 fHistAllRecETDep->Fill(xCluster,etDep);
 279                 fHistAllRec->Fill(xCluster);
 280
 281                 if (track)
 282                 {
 283                         xCharged[0] = track->Eta();
 284                         xCharged[1] = track->Pt();
 285                 }
 286                 else
 287                 {
 288                         xCharged[0] = -99;
 289                         xCharged[1] = -99;
 290                 }
 291                 xCharged[2] = caloE;
 292                 xCharged[3] = caloPos.Eta();
 293                 xCharged[4] = TMath::RadToDeg()*caloPos.Phi();
 294                 xCharged[5] = caloCluster->GetNCells();
 295                 xCharged[6] = res;
 296
 297                 Bool_t isCharged = kFALSE;
 298
 299                 if (maxPid == AliPID::kProton)
 300                 {
 301                         /*
 302                         fHistProtonRecResEDepETDep->Fill(caloE,res,etDep);
 303                         fHistProtonRecResPtETDep->Fill(track->Pt(),res,etDep);
 304                         fHistProtonRecResEDep->Fill(caloE,res,etDep);
 305                         fHistProtonRecResPt->Fill(track->Pt(),res,etDep);
 306                         */
 307
 308                         fHistProtonRecETDep->Fill(xCharged,etDep);
 309                         fHistProtonRec->Fill(xCharged);
 310
 311                         fHistProtonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 312
 313                         if ((res>0.) && (res<fResCut))
 314                         {
 315                                 /*
 316                                 fHistProtonMatchEtaEDepETDep->Fill(caloE,track->Eta(),etDep);
 317                                 fHistProtonMatchEtaPtETDep->Fill(track->Pt(),track->Eta(),etDep);
 318                                 fHistProtonMatchEtaETDep->Fill(etDep,track->Eta());
 319                                 fHistProtonMatchEtaPt->Fill(track->Pt(),track->Eta());
 320                                  */
 321                                 fProtonMatchtotETDep += etDep;
 322
 323                                 isCharged = kTRUE;
 324                         }
 325                 }
 326                 else if (maxPid == AliPID::kPion)
 327                 {
 328                         /*
 329                         fHistPionRecResEDepETDep->Fill(caloE,res,etDep);
 330                         fHistPionRecResPtETDep->Fill(track->Pt(),res,etDep);
 331                         fHistPionRecResEDep->Fill(caloE,res);
 332                         fHistPionRecResPt->Fill(track->Pt(),res);
 333                         */
 334
 335                         fHistPionRecETDep->Fill(xCharged,etDep);
 336                         fHistPionRec->Fill(xCharged);
 337
 338                         fHistPionRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 339
 340                         if ((res>0.) && (res<fResCut))
 341                         {
 342                                 /*
 343                                 fHistPionMatchEtaEDepETDep->Fill(caloE,track->Eta(),etDep);
 344                                 fHistPionMatchEtaPtETDep->Fill(track->Pt(),track->Eta(),etDep);
 345                                 fHistPionMatchEtaETDep->Fill(etDep,track->Eta());
 346                                 fHistPionMatchEtaPt->Fill(track->Pt(),track->Eta());
 347                                 */
 348
 349                                 fPionMatchtotETDep += etDep;
 350                                 isCharged = kTRUE;
 351                         }
 352                 }
 353                 else if (maxPid == AliPID::kKaon)
 354                 {
 355                         /*
 356                         fHistKaonRecResEDepETDep->Fill(caloE,Res,etDep);
 357                         fHistKaonRecResPtETDep->Fill(track->Pt(),res,etDep);
 358                         fHistKaonRecResEDep->Fill(caloE,res);
 359                         fHistKaonRecResPt->Fill(track->Pt(),res);
 360                         */
 361
 362                         fHistKaonRecETDep->Fill(xCharged,etDep);
 363                         fHistKaonRec->Fill(xCharged);
 364
 365                         fHistKaonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 366
 367                         if ((res>0.) && (res<fResCut))
 368                         {
 369                                 /*
 370                                 fHistKaonMatchEtaEDepETDep->Fill(caloE,track->Eta(),etDep);
 371                                 fHistKaonMatchEtaPtETDep->Fill(track->Pt(),track->Eta(),etDep);
 372                                 fHistKaonMatchEtaETDep->Fill(etDep,track->Eta());
 373                                 fHistKaonMatchEtaPt->Fill(track->Pt(),track->Eta());
 374                                 */
 375
 376                                 fKaonMatchtotETDep += etDep;
 377                                 isCharged = kTRUE;
 378                         }
 379                 }
 380                 else if (maxPid == AliPID::kMuon)
 381                 {
 382                         /*
 383                         fHistMuonRecResEDepETDep->Fill(caloE,res,etDep);
 384                         fHistMuonRecResPtETDep->Fill(track->Pt(),res,etDep);
 385                         fHistMuonRecResEDep->Fill(caloE,res);
 386                         fHistMuonRecResPt->Fill(track->Pt(),res);
 387                         */
 388
 389                         fHistMuonRecETDep->Fill(xCharged,etDep);
 390                         fHistMuonRec->Fill(xCharged);
 391
 392                         fHistMuonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 393
 394                         if ((res>0.) && (res<fResCut))
 395                         {
 396                                 /*
 397                                 fHistMuonMatchEtaEDepETDep->Fill(caloE,track->Eta(),etDep);
 398                                 fHistMuonMatchEtaPtETDep->Fill(track->Pt(),track->Eta(),etDep);
 399                                 fHistMuonMatchEtaETDep->Fill(etDep,track->Eta());
 400                                 fHistMuonMatchEtaPt->Fill(track->Pt(),track->Eta());
 401                                 */
 402
 403                                 fMuonMatchtotETDep += etDep;
 404                                 isCharged = kTRUE;
 405                         }
 406                 }
 407                 else if (maxPid == AliPID::kElectron)
 408                 {
 409                         /*
 410                         fHistElectronRec_ResEDep_ETDep->Fill(caloE,res,etDep);
 411                         fHistElectronRec_ResPt_ETDep->Fill(track->Pt(),res,etDep);
 412                         fHistElectronRec_ResEDep->Fill(caloE,res);
 413                         fHistElectronRec_ResPt->Fill(track->Pt(),res);
 414                         */
 415
 416                         fHistElectronRecETDep->Fill(xCharged,etDep);
 417                         fHistElectronRec->Fill(xCharged);
 418
 419                         fHistElectronRecdEdxP->Fill(track->P(),track->GetTPCsignal());
 420
 421                         if ((res>0.) && (res<fResCut))
 422                         {
 423                                 /*
 424                                 fHistElectronMatchEtaEDepETDep->Fill(caloE,track->Eta(),etDep);
 425                                 fHistElectronMatchEtaPtETDep->Fill(track->Pt(),track->Eta(),etDep);
 426                                 fHistElectronMatchEtaETDep->Fill(etDep,track->Eta());
 427                                 fHistElectronMatchEtaPt->Fill(track->Pt(),track->Eta());
 428                                 */
 429
 430                                 fElectronMatchtotETDep += etDep;
 431                                 isCharged = kTRUE;
 432                         }
 433                 }
 434
 435                 if (!isCharged)
 436                 {
 437                         /*
 438                          fHistNeutralRec_EtaE_ET->Fill(caloE,caloPos.Eta(),etDep);
 439                         //fHistNeutralRec_EtaPt_ET->Fill(caloPos.Pt(),caloPos.Eta(),etDep);
 440                         fHistNeutralRec_EtaET->Fill(etDep,caloPos.Eta());
 441                         fHistNeutralRec_EtaE->Fill(caloE,caloPos.Eta());
 442                         //fHistNeutralRec_EtaPt->Fill(caloPos.Pt(),caloPos.Eta());
 443                          */
 444                         fNeutralRectotET += etDep;
 445                 }
 446
 447         } // end of loop over clusters
 448
 449         fTotEMRectotET = fElectronMatchtotETDep + fNeutralRectotET;
 450         fTotChargedMatchtotETDep = fMuonMatchtotETDep + fPionMatchtotETDep + fKaonMatchtotETDep + fProtonMatchtotETDep;
 451         fTotalRectotETDep = fTotEMRectotET + fTotChargedMatchtotETDep;
 452
 453         fHistAllRectotETDep->Fill(fAllRectotETDep);
 454
 455         fHistElectronMatchtotETDep->Fill(fElectronMatchtotETDep);
 456         fHistNeutralRectotET->Fill(fNeutralRectotET);
 457
 458         fHistTotEMRectotET->Fill(fTotEMRectotET);
 459
 460         fHistMuonMatchtotETDep->Fill(fMuonMatchtotETDep);
 461         fHistPionMatchtotETDep->Fill(fPionMatchtotETDep);
 462         fHistKaonMatchtotETDep->Fill(fKaonMatchtotETDep);
 463         fHistProtonMatchtotETDep->Fill(fProtonMatchtotETDep);
 464         fHistTotChargedMatchtotETDep->Fill(fTotChargedMatchtotETDep);
 465
 466         fHistTotalRectotETDep->Fill(fTotalRectotETDep);
 467
 468         delete fGeoUt;
 469         delete caloClusters;
 470
 471         return 0;
 472 }
 473
 474 void AliAnalysisEmEtReconstructed::Init()
 475 { // init
 476     AliAnalysisEt::Init();
 477 }
 478
 479
 480 void AliAnalysisEmEtReconstructed::ResetEventValues()
 481 { // reset event values
 482         AliAnalysisEt::ResetEventValues();
 483
 484         // collision geometry defaults for p+p:
 485         fAllRectotETDep = 0;
 486
 487         fElectronMatchtotETDep = 0;
 488         fNeutralRectotET = 0;
 489
 490         fTotEMRectotET = 0;
 491
 492         fMuonMatchtotETDep = 0; fPionMatchtotETDep = 0; fKaonMatchtotETDep = 0; fProtonMatchtotETDep = 0;
 493         fTotChargedMatchtotETDep = 0;
 494
 495         fTotalRectotETDep = 0;
 496 }
 497
 498
 499 void AliAnalysisEmEtReconstructed::CreateHistograms()
 500 { // histogram related additions
 501         //AliAnalysisEt::CreateHistograms();
 502
 503         //fHistAllRecEtaEDepETDep = CreateEtaEHisto2D("fHistAllRecEtaEDepETDep_","MC E_{T}, all particles","E_{T}(GeV)");
 504         //fHistAllRecEtaETDep = CreateEtaEtHisto2D("fHistAllRecEtaETDep_","MC all particles","#");
 505
 506         fHistAllRecETDep = CreateClusterHistoSparse("fHistAllRecETDep_","E_{T}, all particles");
 507         fHistAllRec = CreateClusterHistoSparse("fHistAllRec_","counts, all particles");
 508         TString histname = "fHistAllRectotETDep_" + fHistogramNameSuffix;
 509         fHistAllRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
 510
 511         /*
 512         fHistElectronMatchEtaEDepETDep = CreateEtaEHisto2D("fHistElectronMatchEtaEDepETDep_","MC E_{T}, primary Electrons, tracking matched","E_{T} dep (GeV)");
 513         fHistElectronMatchEtaPtETDep = CreateEtaPtHisto2D("fHistElectronMatchEtaPtETDep_","MC E_{T}, primary Electrons","E_{T} dep(GeV)");
 514         fHistElectronMatchEtaETDep = CreateEtaEtHisto2D("fHistElectronMatchEtaETDep_","MC primary Electrons","#");
 515         fHistElectronMatchEtaPt = CreateEtaPtHisto2D("fHistElectronMatchEtaPt_","MC E_{T}, primary Electrons","#");
 516          */
 517         fHistElectronRecETDep = CreateChargedPartHistoSparse("fHistElectronRecETDep_","E_{T}, electrons");
 518         fHistElectronRec = CreateChargedPartHistoSparse("fHistElectronRec_","counts, electrons");
 519         histname = "fHistElectronMatchtotETDep_" + fHistogramNameSuffix;
 520         fHistElectronMatchtotETDep = new TH1F(histname.Data(),"total ET, MC primary Electrons",fgNumOfEBins, fgEAxis);
 521         /*
 522         fHistElectronRec_ResEDep_ETDep = CreateResEHisto2D("fHistElectronRec_ResEDep_ETDep_","MC E_{T}, primary Electrons","E_{T} dep (GeV)");
 523         fHistElectronRec_ResPt_ETDep = CreateResPtHisto2D("fHistElectronRec_ResPt_ETDep_","MC E_{T}, primary Electrons","E_{T} dep (GeV)");
 524         fHistElectronRec_ResEDep = CreateResEHisto2D("fHistElectronRec_ResEDep_","MC primary Electrons","#");
 525         fHistElectronRec_ResPt  = CreateResPtHisto2D("fHistElectronRec_ResPt_","MC primary Electrons","#");
 526          */
 527         histname = "fHistElectronRecdEdxP_" + fHistogramNameSuffix;
 528         fHistElectronRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 529
 530         /*
 531         fHistNeutralRec_EtaE_ET = CreateEtaEHisto2D("fHistNeutralRec_EtaE_ET_","MC E_{T}, primary Neutrals","E_{T}(GeV)");
 532         //fHistNeutralRec_EtaPt_ET = CreateEtaPtHisto2D("fHistNeutralRec_EtaPt_ET_","MC E_{T}, primary Neutrals","E_{T}(GeV)");
 533         fHistNeutralRec_EtaET = CreateEtaEtHisto2D("fHistNeutralRec_EtaET_","MC primary Neutrals","#");
 534         fHistNeutralRec_EtaE = CreateEtaEHisto2D("fHistNeutralRec_EtaE_","MC primary Neutrals","#");
 535         //fHistNeutralRec_EtaPt = CreateEtaPtHisto2D("fHistNeutralRec_EtaPt_","MC primary Neutrals","#");
 536         */
 537         histname = "fHistNeutralRectotET_" + fHistogramNameSuffix;
 538         fHistNeutralRectotET = new TH1F(histname.Data(),"total ET, neutral particles",fgNumOfEBins, fgEAxis);
 539
 540         histname = "fHistTotEMRectotET_" + fHistogramNameSuffix;
 541         fHistTotEMRectotET = new TH1F(histname.Data(),"total electromagnetic ET",fgNumOfEBins, fgEAxis);
 542
 543         /*
 544         fHistMuonMatchEtaEDepETDep = CreateEtaEHisto2D("fHistMuonMatchEtaEDepETDep_","MC E_{T}, primary Muons, tracking matched","E_{T} dep (GeV)");
 545         fHistMuonMatchEtaPtETDep = CreateEtaPtHisto2D("fHistMuonMatchEtaPtETDep_","MC E_{T}, primary Muons","E_{T} dep(GeV)");
 546         fHistMuonMatchEtaETDep = CreateEtaEtHisto2D("fHistMuonMatchEtaETDep_","MC primary Muons","#");
 547         fHistMuonMatchEtaPt = CreateEtaPtHisto2D("fHistMuonMatchEtaPt_","MC E_{T}, primary Muons","#");
 548         */
 549         fHistMuonRecETDep = CreateChargedPartHistoSparse("fHistMuonRecETDep_","E_{T}, muons");
 550         fHistMuonRec = CreateChargedPartHistoSparse("fHistMuonRec_","counts, muons");
 551         histname = "fHistMuonMatchtotETDep_" + fHistogramNameSuffix;
 552         fHistMuonMatchtotETDep = new TH1F(histname.Data(),"total ET, Muons",fgNumOfEBins, fgEAxis);
 553         /*
 554         fHistMuonRecResEDepETDep = CreateResEHisto2D("fHistMuonRecResEDepETDep_","MC E_{T}, primary Muons","E_{T} dep (GeV)");
 555         fHistMuonRecResPtETDep = CreateResPtHisto2D("fHistMuonRecResPtETDep_","MC E_{T}, primary Muons","E_{T} dep (GeV)");
 556         fHistMuonRecResEDep = CreateResEHisto2D("fHistMuonRecResEDep_","MC primary Muons","#");
 557         fHistMuonRecResPt  = CreateResPtHisto2D("fHistMuonRecResPt_","MC primary Muons","#");
 558         */
 559         histname = "fHistMuonRecdEdxP_" + fHistogramNameSuffix;
 560         fHistMuonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 561
 562         /*
 563         fHistPionMatchEtaEDepETDep = CreateEtaEHisto2D("fHistPionMatchEtaEDepETDep_","MC E_{T}, primary Pions, tracking matched","E_{T} dep (GeV)");
 564         fHistPionMatchEtaPtETDep = CreateEtaPtHisto2D("fHistPionMatchEtaPtETDep_","MC E_{T}, primary Pions","E_{T} dep(GeV)");
 565         fHistPionMatchEtaETDep = CreateEtaEtHisto2D("fHistPionMatchEtaETDep_","MC primary Pions","#");
 566         fHistPionMatchEtaPt = CreateEtaPtHisto2D("fHistPionMatchEtaPt_","MC E_{T}, primary Pions","#");
 567         */
 568         fHistPionRecETDep = CreateChargedPartHistoSparse("fHistPionRecETDep_","E_{T}, pions");
 569         fHistPionRec = CreateChargedPartHistoSparse("fHistPionRec_","counts, pions");
 570         histname = "fHistPionMatchtotETDep_" + fHistogramNameSuffix;
 571         fHistPionMatchtotETDep = new TH1F(histname.Data(),"total ET, Pions",fgNumOfEBins, fgEAxis);
 572         /*
 573         fHistPionRecResEDepETDep = CreateResEHisto2D("fHistPionRecResEDepETDep_","MC E_{T}, primary Pions","E_{T} dep (GeV)");
 574         fHistPionRecResPtETDep = CreateResPtHisto2D("fHistPionRecResPtETDep_","MC E_{T}, primary Pions","E_{T} dep (GeV)");
 575         fHistPionRecResEDep = CreateResEHisto2D("fHistPionRecResEDep_","MC primary Pions","#");
 576         fHistPionRecResPt  = CreateResPtHisto2D("fHistPionRecResPt_","MC primary Pions","#");
 577         */
 578         histname = "fHistPionRecdEdxP_" + fHistogramNameSuffix;
 579         fHistPionRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 580
 581         /*
 582         fHistKaonMatchEtaEDepETDep = CreateEtaEHisto2D("fHistKaonMatchEtaEDepETDep_","MC E_{T}, primary Kaons, tracking matched","E_{T} dep (GeV)");
 583         fHistKaonMatchEtaPtETDep = CreateEtaPtHisto2D("fHistKaonMatchEtaPtETDep_","MC E_{T}, primary Kaons","E_{T} dep(GeV)");
 584         fHistKaonMatchEtaETDep = CreateEtaEtHisto2D("fHistKaonMatchEtaETDep_","MC primary Kaons","#");
 585         fHistKaonMatchEtaPt = CreateEtaPtHisto2D("fHistKaonMatchEtaPt_","MC primary Kaons","#");
 586         */
 587         fHistKaonRecETDep = CreateChargedPartHistoSparse("fHistKaonRecETDep_","E_{T}, kaons");
 588         fHistKaonRec = CreateChargedPartHistoSparse("fHistKaonRec_","counts, kaons");
 589         histname = "fHistKaonMatchtotETDep_" + fHistogramNameSuffix;
 590         fHistKaonMatchtotETDep = new TH1F(histname.Data(),"total ET, Kaons",fgNumOfEBins, fgEAxis);
 591         /*
 592         fHistKaonRecResEDepETDep = CreateResEHisto2D("fHistKaonRecResEDepETDep_","MC E_{T}, primary Kaons","E_{T} dep (GeV)");
 593         fHistKaonRecResPtETDep = CreateResPtHisto2D("fHistKaonRecResPtETDep_","MC E_{T}, primary Kaons","E_{T} dep (GeV)");
 594         fHistKaonRecResEDep = CreateResEHisto2D("fHistKaonRecResEDep_","MC primary Kaons","#");
 595         fHistKaonRecResPt  = CreateResPtHisto2D("fHistKaonRecResPt_","MC primary Kaons","#");
 596         */
 597         histname = "fHistKaonRecdEdxP_" + fHistogramNameSuffix;
 598         fHistKaonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 599
 600         /*
 601         fHistProtonMatchEtaEDepETDep = CreateEtaEHisto2D("fHistProtonMatchEtaEDepETDep_","MC E_{T}, primary Protons, tracking matched","E_{T} dep (GeV)");
 602         fHistProtonMatchEtaPtETDep = CreateEtaPtHisto2D("fHistProtonMatchEtaPtETDep_","MC E_{T}, primary Protons","E_{T} dep(GeV)");
 603         fHistProtonMatchEtaETDep = CreateEtaEtHisto2D("fHistProtonMatchEtaETDep_","MC primary Protons","#");
 604         fHistProtonMatchEtaPt = CreateEtaPtHisto2D("fHistProtonMatchEtaPt_","MC primary Protons","#");
 605         */
 606         fHistProtonRecETDep = CreateChargedPartHistoSparse("fHistProtonRecETDep_","E_{T}, protons");
 607         fHistProtonRec = CreateChargedPartHistoSparse("fHistProtonRec_","counts, protons");
 608         histname = "fHistProtonMatchtotETDep_" + fHistogramNameSuffix;
 609         fHistProtonMatchtotETDep = new TH1F(histname.Data(),"total ET, Protons",fgNumOfEBins, fgEAxis);
 610         /*
 611         fHistProtonRecResEDepETDep = CreateResEHisto2D("fHistProtonRecResEDepETDep_","MC E_{T}, primary Protons","E_{T} dep (GeV)");
 612         fHistProtonRecResPtETDep = CreateResPtHisto2D("fHistProtonRecResPtETDep_","MC E_{T}, primary Protons","E_{T} dep (GeV)");
 613         fHistProtonRecResEDep = CreateResEHisto2D("fHistProtonRecResEDep_","MC primary Protons","#");
 614         fHistProtonRecResPt  = CreateResPtHisto2D("fHistProtonRecResPt_","MC primary Protons","#");
 615         */
 616         histname = "fHistProtonRecdEdxP_" + fHistogramNameSuffix;
 617         fHistProtonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
 618
 619         histname = "fHistTotChargedMatchtotETDep_" + fHistogramNameSuffix;
 620         fHistTotChargedMatchtotETDep = new TH1F(histname.Data(),"total ET, charged particles",fgNumOfEBins, fgEAxis);
 621
 622         histname = "fHistTotalRectotETDep_" + fHistogramNameSuffix;
 623         fHistTotalRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
 624
 625         histname = "fHistDeltaRZ_" + fHistogramNameSuffix;
 626         fHistDeltaRZ = new TH2F(histname,"#Delta#phi vs #Delta#eta (track projection - cluster position)",200,-0.1,0.1,200,-0.1,0.1);
 627 }
 628
 629 void AliAnalysisEmEtReconstructed::FillOutputList(TList *list)
 630 {//Function for filling the output list
 631         //AliAnalysisEt::FillOutputList(list);
 632
 633         //list->Add(fHistAllRecEtaEDepETDep);
 634         //list->Add(fHistAllRecEtaETDep);
 635         list->Add(fHistAllRecETDep);
 636         list->Add(fHistAllRec);
 637         list->Add(fHistAllRectotETDep);
 638
 639         /*
 640         list->Add(fHistElectronMatchEtaEDepETDep);
 641         list->Add(fHistElectronMatchEtaPtETDep);
 642         list->Add(fHistElectronMatchEtaETDep);
 643         list->Add(fHistElectronMatchEtaPt);
 644
 645         list->Add(fHistElectronRec_ResEDep_ETDep);
 646         list->Add(fHistElectronRec_ResPt_ETDep);
 647         list->Add(fHistElectronRec_ResEDep);
 648         list->Add(fHistElectronRec_ResPt);
 649         */
 650
 651         list->Add(fHistElectronRecETDep);
 652         list->Add(fHistElectronRec);
 653         list->Add(fHistElectronMatchtotETDep);
 654         list->Add(fHistElectronRecdEdxP);
 655
 656         /*
 657         list->Add(fHistNeutralRec_EtaE_ET);
 658         //list->Add(fHistNeutralRec_EtaPt_ET);
 659         list->Add(fHistNeutralRec_EtaET);
 660         list->Add(fHistNeutralRec_EtaE);
 661         //list->Add(fHistNeutralRec_EtaPt);
 662         */
 663         list->Add(fHistNeutralRectotET);
 664
 665         list->Add(fHistTotEMRectotET);
 666
 667         /*
 668         list->Add(fHistMuonMatchEtaEDepETDep);
 669         list->Add(fHistMuonMatchEtaPtETDep);
 670         list->Add(fHistMuonMatchEtaETDep);
 671         list->Add(fHistMuonMatchEtaPt);
 672
 673         list->Add(fHistMuonRecResEDepETDep);
 674         list->Add(fHistMuonRecResPtETDep);
 675         list->Add(fHistMuonRecResEDep);
 676         list->Add(fHistMuonRecResPt);
 677          */
 678         list->Add(fHistMuonRecETDep);
 679         list->Add(fHistMuonRec);
 680         list->Add(fHistMuonRecdEdxP);
 681         list->Add(fHistMuonMatchtotETDep);
 682
 683         /*
 684         list->Add(fHistPionMatchEtaEDepETDep);
 685         list->Add(fHistPionMatchEtaPtETDep);
 686         list->Add(fHistPionMatchEtaETDep);
 687         list->Add(fHistPionMatchEtaPt);
 688
 689         list->Add(fHistPionRecResEDepETDep);
 690         list->Add(fHistPionRecResPtETDep);
 691         list->Add(fHistPionRecResEDep);
 692         list->Add(fHistPionRecResPt);
 693         */
 694         list->Add(fHistPionRecETDep);
 695         list->Add(fHistPionRec);
 696         list->Add(fHistPionMatchtotETDep);
 697         list->Add(fHistPionRecdEdxP);
 698
 699         /*
 700         list->Add(fHistKaonMatchEtaEDepETDep);
 701         list->Add(fHistKaonMatchEtaPtETDep);
 702         list->Add(fHistKaonMatchEtaETDep);
 703         list->Add(fHistKaonMatchEtaPt);
 704
 705         list->Add(fHistKaonRecResEDepETDep);
 706         list->Add(fHistKaonRecResPtETDep);
 707         list->Add(fHistKaonRecResEDep);
 708         list->Add(fHistKaonRecResPt);
 709         */
 710
 711         list->Add(fHistKaonRecETDep);
 712         list->Add(fHistKaonRec);
 713         list->Add(fHistKaonMatchtotETDep);
 714         list->Add(fHistKaonRecdEdxP);
 715
 716         /*
 717         list->Add(fHistProtonMatchEtaEDepETDep);
 718         list->Add(fHistProtonMatchEtaPtETDep);
 719         list->Add(fHistProtonMatchEtaETDep);
 720         list->Add(fHistProtonMatchEtaPt);
 721
 722         list->Add(fHistProtonRecResEDepETDep);
 723         list->Add(fHistProtonRecResPtETDep);
 724         list->Add(fHistProtonRecResEDep);
 725         list->Add(fHistProtonRecResPt);
 726         */
 727
 728         list->Add(fHistProtonRecETDep);
 729         list->Add(fHistProtonRec);
 730         list->Add(fHistProtonMatchtotETDep);
 731         list->Add(fHistProtonRecdEdxP);
 732
 733         list->Add(fHistTotChargedMatchtotETDep);
 734         list->Add(fHistTotalRectotETDep);
 735
 736         list->Add(fHistDeltaRZ);
 737 }
 738
 739 //________________________________________________________________________
 740 // project to a EMCal radius
 741 Bool_t AliAnalysisEmEtReconstructed::GetTrackProjection(AliExternalTrackParam *trackParam, TVector3 &trackPos)
 742 {//Gets the projection of the track
 743     Bool_t proj = kFALSE;
 744     Double_t emcalR = fGeoUt->GetEMCGeometry()->GetIPDistance();
 745
 746     if (trackParam) //it is constructed from TParticle
 747     {
 748         Double_t trkPos[3] = {0};
 749
 750         //Assume the track is a pion with mass 0.139GeV/c^2
 751         //Extrapolation step is 1cm
 752         if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, emcalR, 0.139, 1, kTRUE, 0.8) ) return proj;
 753
 754         trackParam->GetXYZ(trkPos);
 755
 756         trackPos.SetXYZ(trkPos[0],trkPos[1],trkPos[2]);
 757
 758         proj = kTRUE;
 759     }
 760
 761     return proj;
 762 }
 763
 764 //________________________________________________________________________
 765 // project to a cluster position
 766 Bool_t AliAnalysisEmEtReconstructed::GetTrackProjection(AliEMCALTrack* emcTrack, TVector3 &trackPos, TVector3 clusPos)
 767 {//project to a cluster position
 768         Bool_t proj = kFALSE;
 769
 770         if (emcTrack)
 771         {
 772                 Double_t trkPos[3] = {0};
 773
 774                 emcTrack->PropagateToGlobal(clusPos.X(),clusPos.Y(),clusPos.Z(),0.,0.);
 775                 emcTrack->GetXYZ(trkPos);
 776
 777                 trackPos.SetXYZ(trkPos[0],trkPos[1],trkPos[2]);
 778
 779                 proj = kTRUE;
 780         }
 781
 782         return proj;
 783 }
 784
 785 //________________________________________________________________________
 786 AliESDtrack* AliAnalysisEmEtReconstructed::FindMatch(const AliESDCaloCluster *caloCluster, Double_t& resMin)
 787 {//find a matched track
 788         Double_t res=0;
 789         resMin=999;
 790
 791         TVector3 caloPos(0,0,0);
 792         Float_t pos[3] = {0};
 793         caloCluster->GetPosition(pos);
 794         caloPos.SetXYZ(pos[0],pos[1],pos[2]);
 795
 796         // loop over tracks
 797         TVector3 trackPos(0,0,0);
 798         TVector3 trackMatchPos(0,0,0);
 799         AliEMCALTrack *emcTrack = 0;
 800         AliESDtrack *trackMatch = 0;
 801
 802         TObjArray* list = fEsdtrackCutsITSTPC->GetAcceptedTracks(fESD);;
 803         Int_t nGoodTracks = list->GetEntries();
 804
 805         for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
 806         {
 807                 AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
 808         if (!track)
 809         {
 810             AliError(Form("ERROR: Could not get track %d", iTrack));
 811             continue;
 812         }
 813
 814                 emcTrack = new AliEMCALTrack(*track);
 815
 816                 if (GetTrackProjection(emcTrack,trackPos,caloPos))
 817                 {
 818                         res = sqrt(pow(trackPos.Phi()-caloPos.Phi(),2)+pow(trackPos.Eta()-caloPos.Eta(),2));
 819
 820                         if (res < resMin)
 821                         {
 822                                 resMin = res;
 823                                 trackMatch = track;
 824                                 trackMatchPos.SetXYZ(trackPos.X(),trackPos.Y(),trackPos.Z());
 825                         }
 826                 }
 827
 828                 delete emcTrack;
 829         }
 830
 831         fHistDeltaRZ->Fill(trackMatchPos.Phi()-caloPos.Phi(),trackMatchPos.Eta()-caloPos.Eta());
 832
 833         return trackMatch;
 834 }
 835
 836 //________________________________________________________________________
 837 Double_t AliAnalysisEmEtReconstructed::GetTrackPID(const AliESDtrack *track) const
 838 {//Get the default track ID
 839         const Double_t *pidWeights = track->PID();
 840         Int_t maxpid = -1;
 841         Double_t maxpidweight = 0;
 842
 843         if (pidWeights)
 844         {
 845                 for (Int_t p =0; p < AliPID::kSPECIES; p++)
 846                 {
 847                         if (pidWeights[p] > maxpidweight)
 848                         {
 849                                 maxpidweight = pidWeights[p];
 850                                 maxpid = p;
 851                         }
 852                 }
 853         }
 854
 855         return maxpid;
 856 }