PWG4/AliAnaGamma.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16
  17 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.4  2007/10/29 13:48:42  gustavo
  21  * Corrected coding violations
  22  *
  23  * Revision 1.2  2007/08/17 12:40:04  schutz
  24  * New analysis classes by Gustavo Conesa
  25  *
  26  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  27  * new analysis classes in the the new analysis framework
  28  *
  29  *
  30  */
  31
  32 //_________________________________________________________________________
  33 // Base class for gamma and correlation analysis
  34 // It is called by the task class AliAnalysisGammaTask and it connects the input (ESD/AOD/MonteCarlo)
  35 // got with AliGammaReader (produces TClonesArrays of TParticles), with the analysis classes
  36 // AliAnaGammaDirect, AliAnaGammaCorrelation ....
  37 //
  38 //*-- Author: Gustavo Conesa (INFN-LNF)
  39
  40 // --- ROOT system ---
  41
  42 #include <TParticle.h>
  43 #include <TH2.h>
  44
  45 //---- AliRoot system ----
  46 #include "AliAnaGamma.h"
  47 #include "AliGammaReader.h"
  48 #include "AliAnaGammaDirect.h"
  49 #include "AliAnaGammaCorrelation.h"
  50 #include "AliAnaGammaSelection.h"
  51 #include "AliNeutralMesonSelection.h"
  52 #include "AliAODCaloCluster.h"
  53 #include "AliAODEvent.h"
  54 #include "Riostream.h"
  55 #include "AliLog.h"
  56
  57 ClassImp(AliAnaGamma)
  58
  59
  60 //____________________________________________________________________________
  61   AliAnaGamma::AliAnaGamma() :
  62     TObject(),
  63     fOutputContainer(0x0),
  64     fAnaType(0),  fCalorimeter(0), fData(0x0), fKine(0x0),
  65     fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0), fGammaSelection(0x0),
  66     fNeutralMesonSelection(0x0), fAODclusters(0x0), fNAODclusters(0)
  67 {
  68   //Default Ctor
  69
  70   //Initialize parameters, pointers and histograms
  71   if(!fReader)
  72     fReader = new AliGammaReader();
  73   if(!fGammaDirect)
  74     fGammaDirect = new AliAnaGammaDirect();
  75   if(!fGammaCorrelation)
  76     fGammaCorrelation = new AliAnaGammaCorrelation();
  77   if(!fGammaSelection)
  78     fGammaSelection = new AliAnaGammaSelection();
  79   if(!fNeutralMesonSelection)
  80     fNeutralMesonSelection = new AliNeutralMesonSelection();
  81
  82   fAODclusters = 0;
  83
  84   InitParameters();
  85
  86 }
  87
  88 //____________________________________________________________________________
  89 AliAnaGamma::AliAnaGamma(const AliAnaGamma & g) :
  90   TObject(),
  91   fOutputContainer(g. fOutputContainer),
  92   fAnaType(g.fAnaType),  fCalorimeter(g.fCalorimeter),
  93   fData(g.fData), fKine(g.fKine),fReader(g.fReader),
  94   fGammaDirect(g.fGammaDirect), fGammaCorrelation(g.fGammaCorrelation),
  95   fGammaSelection(g.fGammaSelection),
  96   fNeutralMesonSelection(g.fNeutralMesonSelection),
  97   fAODclusters(g. fAODclusters), fNAODclusters(g.fNAODclusters)
  98 {
  99   // cpy ctor
 100
 101 }
 102
 103 //_________________________________________________________________________
 104 AliAnaGamma & AliAnaGamma::operator = (const AliAnaGamma & source)
 105 {
 106   // assignment operator
 107
 108   if(this == &source)return *this;
 109   ((TObject *)this)->operator=(source);
 110
 111   fOutputContainer = source.fOutputContainer ;
 112   fAnaType = source.fAnaType;
 113   fCalorimeter = source.fCalorimeter ;
 114   fData = source.fData ;
 115   fKine = source.fKine ;
 116   fReader = source.fReader ;
 117   fGammaDirect = source.fGammaDirect ;
 118   fGammaCorrelation = source.fGammaCorrelation ;
 119   fGammaSelection = source.fGammaSelection ;
 120   fNeutralMesonSelection = source.fNeutralMesonSelection ;
 121
 122   return *this;
 123
 124 }
 125
 126 //____________________________________________________________________________
 127 AliAnaGamma::~AliAnaGamma()
 128 {
 129   // Remove all pointers.
 130
 131   // Protection added in case of NULL pointers (MG)
 132   if (fOutputContainer) {
 133      fOutputContainer->Clear();
 134      delete fOutputContainer ;
 135   }
 136
 137   if (fData) delete fData ;
 138   if (fKine) delete fKine ;
 139   if (fReader) delete fReader ;
 140   if (fGammaDirect) delete fGammaDirect ;
 141   if (fGammaCorrelation) delete fGammaCorrelation ;
 142   if (fGammaSelection) delete fGammaSelection ;
 143   if (fNeutralMesonSelection) delete fNeutralMesonSelection ;
 144
 145 }
 146
 147 //________________________________________________________________________
 148 void AliAnaGamma::Init()
 149 {
 150
 151   //Init container histograms and other common variables
 152
 153   //Histograms container
 154   fOutputContainer = new TList ;
 155
 156   //Fill container with appropriate histograms
 157
 158   //Set selection  analysis histograms
 159   TList * selectcontainer =  fGammaSelection->GetCreateOutputObjects();
 160   for(Int_t i = 0; i < selectcontainer->GetEntries(); i++){
 161     Bool_t  add = kTRUE ;
 162     TString name = (selectcontainer->At(i))->GetName();
 163     if(!fReader->IsEMCALOn() && name.Contains("EMCAL")) add = kFALSE;
 164     if(!fReader->IsPHOSOn() && name.Contains("PHOS"))   add = kFALSE;
 165     if(!fReader->IsCTSOn() &&  !fGammaSelection->FillCTS() && name.Contains("CTS"))   add = kFALSE;
 166     if(add) fOutputContainer->Add(selectcontainer->At(i)) ;
 167   }  //Set selection  analysis histograms
 168
 169
 170   //Set prompt photon analysis histograms
 171   TList * promptcontainer =  fGammaDirect->GetCreateOutputObjects();
 172   for(Int_t i = 0; i < promptcontainer->GetEntries(); i++)
 173     fOutputContainer->Add(promptcontainer->At(i)) ;
 174
 175   //Check if selected options are correct or set them when necessary
 176   if(fReader->GetDataType() == AliGammaReader::kMCData){
 177     fGammaDirect->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
 178     fGammaSelection->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
 179   }
 180
 181   if(fAnaType == kCorrelation){
 182
 183     //Check if selected options are correct
 184     if (fGammaDirect->GetICMethod()==AliAnaGammaDirect::kSeveralIC)
 185       AliFatal("Correlation not allowed with multiple isolation cuts, kCorrelation and kSeveralIC do not go together");
 186
 187     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kParton &&
 188        fReader->GetDataType() != AliGammaReader::kMC)
 189       AliFatal("kParton must be analyzed with data kMC");
 190
 191     //Set the parameters for the neutral pair selection depending on the analysis,
 192     fNeutralMesonSelection->SetDeltaPhiCutRange(fGammaCorrelation->GetDeltaPhiMinCut(),
 193                                                 fGammaCorrelation->GetDeltaPhiMaxCut());
 194
 195
 196     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kHadron){
 197       fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiMinPt);
 198       fNeutralMesonSelection->SetMinPt(fGammaCorrelation->GetMinPtHadron());
 199
 200     }
 201
 202     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kJetLeadCone){
 203       fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiPtRatio);
 204       fNeutralMesonSelection->SetRatioCutRange(fGammaCorrelation->GetRatioMinCut(),
 205                                                fGammaCorrelation->GetRatioMaxCut());
 206     }
 207
 208     //Set the neutral mesosn selection histograms
 209     TList * neutralmesoncontainer =  fNeutralMesonSelection->GetCreateOutputObjects();
 210     if(fNeutralMesonSelection->AreNeutralMesonSelectionHistosKept()){
 211       for(Int_t i = 0; i < neutralmesoncontainer->GetEntries(); i++)
 212         fOutputContainer->Add(neutralmesoncontainer->At(i)) ;
 213     }
 214
 215     //Set correlation histograms
 216     TList * correlationcontainer =  fGammaCorrelation->GetCreateOutputObjects();
 217     for(Int_t i = 0; i < correlationcontainer->GetEntries(); i++)
 218       fOutputContainer->Add(correlationcontainer->At(i)) ;
 219     fGammaCorrelation->SetOutputContainer(fOutputContainer);
 220     fGammaCorrelation->SetNeutralMesonSelection(fNeutralMesonSelection);
 221
 222   }//kCorrelation
 223
 224 }
 225
 226 //____________________________________________________________________________
 227 void AliAnaGamma::InitParameters()
 228 {
 229
 230   //Init data members
 231
 232   fAnaType = kPrompt;
 233   fCalorimeter = "EMCAL";
 234
 235 }
 236
 237 //__________________________________________________________________
 238 void AliAnaGamma::Print(const Option_t * opt) const
 239 {
 240
 241   //Print some relevant parameters set for the analysis
 242   if(! opt)
 243     return;
 244
 245   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 246   printf("Analysis type           =     %d\n", fAnaType) ;
 247   printf("Calorimeter           =     %s\n", fCalorimeter.Data()) ;
 248
 249   switch(fAnaType)
 250     {
 251     case kPrompt:
 252       {
 253         fGammaDirect->Print("");
 254       }// case kIsolationCut
 255       break;
 256
 257     case kCorrelation:
 258       {
 259         fGammaCorrelation->Print("");
 260       }//  case kCorrelation
 261       break;
 262
 263     }//switch
 264 }
 265
 266
 267 //____________________________________________________________________________
 268 Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){
 269
 270   AliDebug(1,Form("Entry %d",entry));
 271   //cout<<"Event >>>>>>>>>>>>> "<<entry<<endl;
 272   if(!fOutputContainer)
 273     AliFatal("Histograms not initialized");
 274
 275   //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
 276
 277   TClonesArray * plCTS      = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
 278   TClonesArray * plEMCAL    = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
 279   TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured  Gamma calorimeter
 280   TClonesArray * plPrimCTS      = new TClonesArray("TParticle",1000); // primary tracks
 281   TClonesArray * plPrimEMCAL    = new TClonesArray("TParticle",1000);   // primary emcal clusters
 282   TClonesArray * plPrimPHOS     = new TClonesArray("TParticle",1000);  // primary phos clusters
 283   TClonesArray * plParton   = new TClonesArray("TParticle",1000);  // All partons
 284   //Fill lists with photons, neutral particles and charged particles
 285   //look for the highest energy photon in the event inside fCalorimeter
 286
 287   //Fill particle lists
 288   if(fReader->GetDataType() == AliGammaReader::kData){
 289     AliDebug(1,"Data analysis");
 290     fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL,NULL,NULL);
 291   }
 292   else if( fReader->GetDataType()== AliGammaReader::kMC){
 293     AliDebug(1,"Kinematics analysis");
 294     fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton,NULL,NULL);
 295   }
 296   else if(fReader->GetDataType() == AliGammaReader::kMCData) {
 297    AliDebug(1,"Data + Kinematics analysis");
 298    fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,plPrimCTS,plPrimEMCAL,plPrimPHOS);
 299   }
 300   else
 301     AliError("Option not implemented");
 302
 303   //Fill AOD with calorimeter
 304   //Temporal solution, just for testing
 305   //FillAODs(plPHOS,plEMCAL);
 306
 307   //Select particles to do the final analysis.
 308   if(fReader->IsEMCALOn()) fGammaSelection->Selection("EMCAL",plEMCAL,plPrimEMCAL);
 309   if(fReader->IsPHOSOn())  fGammaSelection->Selection("PHOS",plPHOS,plPrimPHOS);
 310   if(fReader->IsCTSOn() &&  fGammaSelection->FillCTS()) fGammaSelection->Selection("CTS",plCTS,plPrimCTS);
 311
 312
 313   //Search highest energy prompt gamma in calorimeter
 314   if(fCalorimeter == "PHOS")
 315     MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ;
 316   else if (fCalorimeter == "EMCAL")
 317     MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton, plPrimEMCAL) ;
 318   else
 319     AliFatal("Wrong calorimeter name");
 320
 321   plCTS->Clear() ;
 322   plEMCAL->Clear() ;
 323   plPHOS->Clear() ;
 324   plParton->Clear() ;
 325   plPrimCTS->Clear() ;
 326   plPrimEMCAL->Clear() ;
 327   plPrimPHOS->Clear() ;
 328
 329   delete plCTS ;
 330   delete plPHOS ;
 331   delete plEMCAL ;
 332   delete plParton ;
 333   delete plPrimCTS ;
 334   delete plPrimPHOS ;
 335   delete plPrimEMCAL ;
 336
 337   return kTRUE;
 338
 339 }
 340
 341 //____________________________________________________________________________
 342 void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton, TClonesArray * plPrimCalo)  {
 343
 344   TParticle * pGamma = new TParticle ;
 345   Bool_t isInCalo = kFALSE ;
 346
 347   switch(fAnaType)
 348     {
 349
 350       //Only Prompt photon analysis
 351     case kPrompt:
 352       {
 353         AliDebug(1,"kPrompt analysis");
 354         switch(fGammaDirect->GetICMethod())
 355           {
 356
 357           case AliAnaGammaDirect::kSeveralIC:
 358             {
 359               fGammaDirect->MakeSeveralICAnalysis(plCalo, plCTS);
 360               AliDebug(1,"kSeveralIC analysis");
 361             }
 362             break;
 363
 364           default :
 365             {
 366               fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
 367               if(!isInCalo)
 368                 AliDebug(1,"Prompt gamma not found");
 369             }
 370             break;
 371           }//IC method
 372       }// case kPrompt:
 373       break;
 374
 375       //Correlate prompt photon with something: parton, hadron, jet.
 376     case kCorrelation:
 377       {
 378         AliDebug(1,"kCorrelation analysis");
 379         //Find prompt photon
 380         fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
 381
 382         if(isInCalo){//If prompt photon found, do correlation
 383
 384           switch(fGammaCorrelation->GetCorrelationType())
 385             {
 386             case AliAnaGammaCorrelation::kParton:
 387               {
 388                 AliDebug(1,"kParton correlation");
 389                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plParton, NULL);
 390               }//  case kParton
 391               break;
 392
 393             case AliAnaGammaCorrelation::kHadron:
 394               {
 395                 AliDebug(1,"kHadron correlation");
 396                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
 397               }//  case kHadron
 398               break;
 399
 400             case AliAnaGammaCorrelation::kJetLeadCone:
 401               {
 402                 AliDebug(1,"kJetLeadCone correlation");
 403                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
 404               }//  case kJetLeadCone
 405               break;
 406
 407             case AliAnaGammaCorrelation::kJetFinder:
 408               {
 409                 AliDebug(1,"kJetFinder correlation");
 410                 printf("Analysis not implemented \n");
 411               }//  case kJetFinder
 412               break;
 413             }// switch correlation
 414         }// is in calo
 415         else  AliDebug(2,"Prompt gamma not found");
 416       }// case kCorrelation
 417       break;
 418
 419     } //switch(fAnaType)
 420
 421   delete pGamma ;
 422
 423 }
 424
 425 //____________________________________________________
 426 void AliAnaGamma::AddCluster(AliAODCaloCluster p)
 427 {
 428 // Add new jet to the list
 429   cout<<"AOD list pointer "<<fAODclusters<<" nAODclusters "<<fNAODclusters<<endl;
 430   cout<<"list entries "<<fAODclusters->GetEntries()<<endl;
 431   new ((*fAODclusters)[fNAODclusters++]) AliAODCaloCluster(p);
 432 }
 433
 434 //___________________________________________________
 435 void AliAnaGamma::ConnectAOD(AliAODEvent* aod)
 436 {
 437 // Connect to the AOD
 438   fAODclusters = aod->GetCaloClusters();
 439 }
 440
 441 //____________________________________________________________________________
 442 void AliAnaGamma::FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL){
 443   //Fill AOD caloClusters
 444   //Temporal method, just for testing AOD creation
 445
 446   //Fill AOD with PHOS clusters
 447   Int_t nphos =  plPHOS->GetEntries() ;
 448   cout<<"PHOS entries "<<nphos<<endl;
 449   for(Int_t ipr = 0;ipr < nphos ; ipr ++ ){
 450     TParticle * particle = dynamic_cast<TParticle *>(plPHOS->At(ipr)) ;
 451     Float_t pos []= {0,0,0};
 452     AliAODCaloCluster phos(ipr,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kPHOSNeutral,0);
 453     AddCluster(phos);
 454   }
 455
 456   //Fill AOD with EMCAL clusters
 457   cout<<"EMCAL entries "<<plEMCAL->GetEntries()<<endl;
 458   for(Int_t ipr = 0;ipr < plEMCAL->GetEntries() ; ipr ++ ){
 459     TParticle * particle = dynamic_cast<TParticle *>(plEMCAL->At(ipr)) ;
 460     Float_t pos []= {0,0,0};
 461     AliAODCaloCluster emcal(ipr+nphos,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kEMCALClusterv1,0);
 462     AddCluster(emcal);
 463   }
 464 }