PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalTriggerInfoQA.cxx

   1 #include "AliAnalysisTaskEmcalTriggerInfoQA.h"
   2
   3 #include <Riostream.h>
   4 #include <ctime>
   5 #include <TString.h>
   6 #include <TChain.h>
   7 #include <TTree.h>
   8 #include <TH1D.h>
   9 #include <TH2D.h>
  10 #include <TH3D.h>
  11 #include <TCanvas.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TProfile.h>
  15 #include <TProfile2D.h>
  16 #include <TProfile3D.h>
  17 #include <TRandom.h>
  18 #include <TRandom3.h>
  19
  20 #include "AliAnalysisTaskEmcal.h"
  21 #include "AliAnalysisManager.h"
  22 #include "AliStack.h"
  23 #include "AliESDtrackCuts.h"
  24 #include "AliESDEvent.h"
  25 #include "AliESDInputHandler.h"
  26 #include "AliAODEvent.h"
  27 #include "AliMCEvent.h"
  28 #include "AliEmcalJet.h"
  29 #include "AliEMCALGeometry.h"
  30
  31 #include "AliEmcalTriggerPatchInfo.h"
  32 #include "AliEmcalTriggerSetupInfo.h"
  33
  34 ClassImp(AliAnalysisTaskEmcalTriggerInfoQA)
  35
  36 //________________________________________________________________________
  37 AliAnalysisTaskEmcalTriggerInfoQA::AliAnalysisTaskEmcalTriggerInfoQA() :
  38     AliAnalysisTaskEmcal(),
  39     fOutput(0),
  40     fHistos(0),
  41     fTriggersInfo(0),
  42     fTriggerSetup(0),
  43     fIsInitialized(kFALSE),
  44     fCaloTriggerPatchInfoName("EmcalTriggers"),
  45     fCaloTriggerSetupInfoName("EmcalTriggerSetup")
  46 {
  47 }
  48
  49 //________________________________________________________________________
  50 AliAnalysisTaskEmcalTriggerInfoQA::AliAnalysisTaskEmcalTriggerInfoQA(const char *name) :
  51     AliAnalysisTaskEmcal(name),
  52     fOutput(0),
  53     fHistos(0),
  54     fTriggersInfo(0),
  55     fTriggerSetup(0),
  56     fIsInitialized(kFALSE),
  57     fCaloTriggerPatchInfoName("EmcalTriggers"),
  58     fCaloTriggerSetupInfoName("EmcalTriggerSetup")
  59 {
  60     DefineOutput(1,TList::Class());  // for output list
  61 }
  62
  63 //________________________________________________________________________
  64 AliAnalysisTaskEmcalTriggerInfoQA::~AliAnalysisTaskEmcalTriggerInfoQA()
  65 {
  66 //   // Destructor. Clean-up the output list, but not the histograms that are put inside
  67 //   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
  68 //     if (fHistos && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
  69 //     {
  70 //         delete fOutput;
  71 //     }
  72 //     delete fTrackCuts;
  73 }
  74
  75 //________________________________________________________________________
  76 void AliAnalysisTaskEmcalTriggerInfoQA::UserCreateOutputObjects()
  77 {
  78   Int_t i;
  79
  80   // Create histograms
  81   // Called once (on the worker node)
  82   fIsInitialized=kFALSE;
  83   fOutput = new TList();
  84   fOutput->SetOwner();  // IMPORTANT!
  85
  86   char buf[100];
  87   char strId[3][10];
  88   sprintf( strId[0], "JetLow" );
  89   sprintf( strId[1], "JetHigh" );
  90   sprintf( strId[2], "JetMain" );
  91
  92   fHistos = new TH1*[5000];
  93
  94   for( i = 0; i < 100; i++ )
  95     fHistos[i] = 0;
  96
  97   sprintf( buf, "fhPatchNTotal" );
  98   fHistos[0] = new TH1D( buf, buf, 200, -0.5, 199.5 );
  99   fHistos[0]->SetXTitle( "N_{patch} in event" );
 100   fHistos[0]->SetYTitle( "N" );
 101
 102   // common histos
 103   for( i = 0; i < 3; i++ ){
 104
 105     sprintf( buf, "fhPatchN%s", strId[i] );
 106     fHistos[1000*(i+1)+0] = new TH1D( buf, buf, 200, -0.5, 199.5 );
 107     fHistos[1000*(i+1)+0]->SetXTitle( "N_{patch} in event" );
 108     fHistos[1000*(i+1)+0]->SetYTitle( "N" );
 109
 110     sprintf( buf, "fhThreshold%s", strId[i] );
 111     fHistos[1000*(i+1)+1] = new TH1D( buf, buf, 500, -0.5, 2999.5 );
 112     fHistos[1000*(i+1)+1]->SetXTitle( "thresholds [ADC]" );
 113     fHistos[1000*(i+1)+1]->SetYTitle( "N" );
 114
 115     sprintf( buf, "fhThresholdGeV%s", strId[i] );
 116     fHistos[1000*(i+1)+2] = new TH1D( buf, buf, 500, 0, 250 );
 117     fHistos[1000*(i+1)+2]->SetXTitle( "thresholds [GeV]" );
 118     fHistos[1000*(i+1)+2]->SetYTitle( "N" );
 119
 120     sprintf( buf, "fhADCPatch%s", strId[i] );
 121     fHistos[1000*(i+1)+3] = new TH1D( buf, buf, 500, -0.5, 2999.5 );
 122     fHistos[1000*(i+1)+3]->SetXTitle( "E [ADC]" );
 123     fHistos[1000*(i+1)+3]->SetYTitle( "N" );
 124
 125     sprintf( buf, "fhADCPatchGeV%s", strId[i] );
 126     fHistos[1000*(i+1)+4] = new TH1D( buf, buf, 500, 0, 250 );
 127     fHistos[1000*(i+1)+4]->SetXTitle( "E [GeV]" );
 128     fHistos[1000*(i+1)+4]->SetYTitle( "N" );
 129
 130     sprintf( buf, "fhECells%s", strId[i] );
 131     fHistos[1000*(i+1)+5] = new TH1D( buf, buf, 500, 0, 250 );
 132     fHistos[1000*(i+1)+5]->SetXTitle( "E cells [GeV]" );
 133     fHistos[1000*(i+1)+5]->SetYTitle( "N" );
 134
 135     sprintf( buf, "fhPatchSpacial%s", strId[i] );
 136     fHistos[1000*(i+1)+6] = new TH2D( buf, buf, 9, -0.5, 8.5,
 137                                                 13, -0.5, 12.5 );
 138     fHistos[1000*(i+1)+6]->SetXTitle( "cells in -eta [8 cells]" );
 139     fHistos[1000*(i+1)+6]->SetYTitle( "cells in phi [8 cells]" );
 140     fHistos[1000*(i+1)+6]->SetOption( "COLZ" );
 141     fHistos[1000*(i+1)+6]->SetStats( kFALSE );
 142
 143     sprintf( buf, "fhPatchSpacialADC%s", strId[i] );
 144     fHistos[1000*(i+1)+7] = new TH2D( buf, buf, 9, -0.5, 8.5,
 145                                                 13, -0.5, 12.5 );
 146     fHistos[1000*(i+1)+7]->SetXTitle( "cells in -eta [8 cells]" );
 147     fHistos[1000*(i+1)+7]->SetYTitle( "cells in phi [8 cells]" );
 148     fHistos[1000*(i+1)+7]->SetOption( "COLZ" );
 149     fHistos[1000*(i+1)+7]->SetStats( kFALSE );
 150
 151     sprintf( buf, "fhPatchSpacialECells%s", strId[i] );
 152     fHistos[1000*(i+1)+8] = new TH2D( buf, buf, 9, -0.5, 8.5,
 153                                                 13, -0.5, 12.5 );
 154     fHistos[1000*(i+1)+8]->SetXTitle( "cells in -eta [8 cells]" );
 155     fHistos[1000*(i+1)+8]->SetYTitle( "cells in phi [8 cells]" );
 156     fHistos[1000*(i+1)+8]->SetOption( "COLZ" );
 157     fHistos[1000*(i+1)+8]->SetStats( kFALSE );
 158
 159     sprintf( buf, "fhPatchCenter%s", strId[i] );
 160     fHistos[1000*(i+1)+9] = new TH2D( buf, buf, 100, -0.7, 0.7,
 161                                                 100, 0, TMath::Pi()*2 );
 162     fHistos[1000*(i+1)+9]->SetXTitle( "#eta" );
 163     fHistos[1000*(i+1)+9]->SetYTitle( "#phi [rad]" );
 164     fHistos[1000*(i+1)+9]->SetOption( "COLZ" );
 165     fHistos[1000*(i+1)+9]->SetStats( kFALSE );
 166
 167     sprintf( buf, "fhPatchCenterMass%s", strId[i] );
 168     fHistos[1000*(i+1)+10] = new TH2D( buf, buf, 100, -0.7, 0.7,
 169                                                 100, 0, TMath::Pi()*2 );
 170     fHistos[1000*(i+1)+10]->SetXTitle( "#eta" );
 171     fHistos[1000*(i+1)+10]->SetYTitle( "#phi [rad]" );
 172     fHistos[1000*(i+1)+10]->SetOption( "COLZ" );
 173     fHistos[1000*(i+1)+10]->SetStats( kFALSE );
 174
 175     sprintf( buf, "fhEInPatchSpacial%s", strId[i] );
 176     fHistos[1000*(i+1)+11] = new TH2D( buf, buf, 32, -0.5, 31.5,
 177                                                 32, -0.5, 31.5 );
 178     fHistos[1000*(i+1)+11]->SetXTitle( "cells in -eta" );
 179     fHistos[1000*(i+1)+11]->SetYTitle( "cells in phi" );
 180     fHistos[1000*(i+1)+11]->SetOption( "COLZ" );
 181     fHistos[1000*(i+1)+11]->SetStats( kFALSE );
 182
 183     sprintf( buf, "fhADCvsEcells%s", strId[i] );
 184     fHistos[1000*(i+1)+12] = new TH2D( buf, buf,  500, -0.5, 2999.5,
 185                                                 500, 0, 250 );
 186     fHistos[1000*(i+1)+12]->SetXTitle( "E [ADC]" );
 187     fHistos[1000*(i+1)+12]->SetYTitle( "E cells [GeV]" );
 188     fHistos[1000*(i+1)+12]->SetOption( "COLZ" );
 189     fHistos[1000*(i+1)+12]->SetStats( kFALSE );
 190
 191   }
 192
 193   for( i = 0; i < 5000; i++ )
 194     if( fHistos[i] != 0 )
 195       fOutput->Add( fHistos[i] );
 196
 197
 198   // Post data for ALL output slots >0 here,
 199   // To get at least an empty histogram
 200   // 1 is the outputnumber of a certain weg of task 1
 201   PostData(1, fOutput);
 202 }
 203
 204 void AliAnalysisTaskEmcalTriggerInfoQA::UserExecOnce()
 205 {
 206   AliAnalysisTaskEmcal::ExecOnce();
 207
 208   // Get the event tracks from PicoTracks
 209   fTriggersInfo =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fCaloTriggerPatchInfoName.Data() ));
 210
 211   fTriggerSetup = dynamic_cast <AliEmcalTriggerSetupInfo*>(InputEvent()->FindListObject( fCaloTriggerSetupInfoName.Data() ));
 212
 213   fIsInitialized=kTRUE;
 214 }
 215 //________________________________________________________________________
 216 void AliAnalysisTaskEmcalTriggerInfoQA::UserExec(Option_t *)
 217 {
 218   Int_t nPatch, iPatch, iMainPatch, nJetLow, nJetHigh, type;
 219   AliEmcalTriggerPatchInfo *patch;
 220
 221   if (fIsInitialized==kFALSE)
 222     UserExecOnce();
 223
 224
 225   if( !fTriggersInfo ){
 226     AliError(Form("Calo triggers info container %s not available.", fCaloTriggerPatchInfoName.Data()));
 227     return;
 228   }
 229   if( !fTriggerSetup ){
 230     AliError(Form("Calo trigger setup info container %s not available.", fCaloTriggerSetupInfoName.Data()));
 231     return;
 232   }
 233   if( !fCaloCells ){
 234     AliError(Form("Calo cells container %s not available.", fCaloCellsName.Data()));
 235     return;
 236   }
 237
 238   // first get the thresholds
 239   fHistos[1000+1]->Fill( fTriggerSetup->GetThresholdJetLow() );
 240   fHistos[2000+1]->Fill( fTriggerSetup->GetThresholdJetHigh() );
 241   fHistos[1000+2]->Fill( fTriggerSetup->GetThresholdGeVRoughJetLow() );
 242   fHistos[2000+2]->Fill( fTriggerSetup->GetThresholdGeVRoughJetHigh() );
 243
 244   // now go through patchs
 245   nPatch = fTriggersInfo->GetEntries();
 246
 247   fHistos[0]->Fill( nPatch );
 248
 249   nJetLow = 0;
 250   nJetHigh = 0;
 251   iMainPatch = -1;
 252
 253   for( iPatch = 0; iPatch < nPatch; iPatch++ ){
 254
 255     patch = (AliEmcalTriggerPatchInfo*)fTriggersInfo->At( iPatch );
 256
 257     // check if high/low threshold
 258     // high overrides the low, to avoid double counting
 259     type = -1;
 260     if( patch->IsJetHigh() ){
 261       type = 1;
 262       nJetHigh++;
 263     }
 264     else if( patch->IsJetLow() ){
 265       type = 0;
 266       nJetLow++;
 267     }
 268
 269     if( type == -1 )
 270       continue;
 271
 272     FillPatch( patch, type );
 273
 274     // save main pach position
 275     if( patch->IsMainTrigger() )
 276       iMainPatch = iPatch;
 277
 278   } // patches
 279
 280   // fill jet counts per event
 281   fHistos[1000]->Fill( nJetLow );
 282   fHistos[2000]->Fill( nJetHigh );
 283
 284   //  main trigger patch check -----------------------------
 285   if( iMainPatch == -1 ){
 286       // count of patches
 287       fHistos[3000]->Fill( 0 );
 288   }
 289   else{
 290     patch = (AliEmcalTriggerPatchInfo*)fTriggersInfo->At( iMainPatch );
 291
 292     fHistos[3000]->Fill( 1 );
 293
 294     // what threshold was the main trigger taken with
 295     if( patch->IsJetHigh() ){
 296       fHistos[3000+1]->Fill( fTriggerSetup->GetThresholdJetHigh() );
 297       fHistos[3000+2]->Fill( fTriggerSetup->GetThresholdGeVRoughJetHigh() );
 298     }
 299     else{
 300       fHistos[3000+1]->Fill( fTriggerSetup->GetThresholdJetLow() );
 301       fHistos[3000+2]->Fill( fTriggerSetup->GetThresholdGeVRoughJetLow() );
 302     }
 303
 304     FillPatch( patch, 2 );
 305   }
 306
 307   PostData(1, fOutput);
 308 }
 309
 310
 311 //________________________________________________________________________
 312 void AliAnalysisTaskEmcalTriggerInfoQA::FillPatch( AliEmcalTriggerPatchInfo *patch, Int_t type ){
 313
 314   // fills the patch parameters
 315
 316   Int_t globCol, globRow, i, j, k, absId, cellAbsId[4];;
 317
 318   // patch energies
 319   fHistos[1000*(type+1)+3]->Fill( patch->GetADCAmp() );
 320   fHistos[1000*(type+1)+4]->Fill( patch->GetADCAmpGeVRough() );
 321   fHistos[1000*(type+1)+5]->Fill( patch->GetPatchE() );
 322
 323   //cout << "amp: " << patch->GetADCAmp() << endl;
 324
 325   // get corner, convert from cells to trigger channels
 326   globCol = patch->GetEdgeCellX() / 2;
 327   globRow = patch->GetEdgeCellY() / 2;
 328
 329   // fill in patch steps (8 cells = 4 trigger channels)
 330   ((TH2D*)fHistos[1000*(type+1)+6])->Fill( globCol/4, globRow/4 );
 331   ((TH2D*)fHistos[1000*(type+1)+7])->Fill( globCol/4, globRow/4, patch->GetADCAmp() );
 332   ((TH2D*)fHistos[1000*(type+1)+8])->Fill( globCol/4, globRow/4, patch->GetPatchE() );
 333
 334   ((TH2D*)fHistos[1000*(type+1)+12])->Fill( patch->GetADCAmp(), patch->GetPatchE() );
 335
 336   // phi/eta
 337   ((TH2D*)fHistos[1000*(type+1)+9])->Fill( patch->GetEtaGeo(), patch->GetPhiGeo() );
 338   ((TH2D*)fHistos[1000*(type+1)+10])->Fill( patch->GetEtaCM(), patch->GetPhiCM() );
 339
 340   // E inside patch distribution
 341   // get the absolute trigger ID
 342   fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol, globRow, absId );
 343   // convert to the 4 absId of the cells composing the trigger channel
 344   fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 345
 346   // sum the available energy in the 32/32 window of cells
 347   // step over trigger channels and get all the corresponding cells
 348   for( i = 0; i < 16; i++ ){
 349     for( j = 0; j < 16; j++ ){
 350       // get the 4 cells composing the trigger channel
 351       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+i, globRow+j, absId );
 352       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 353       // add amplitudes and find patch edges
 354       for( k = 0; k < 4; k++ ){
 355         ((TH2D*)fHistos[1000*(type+1)+11])->Fill( i*2+k%2, j*2+k/2,
 356                           fCaloCells->GetCellAmplitude( cellAbsId[k] ));
 357       }
 358     }
 359   } // 32x32 cell window
 360
 361 }
 362
 363 //________________________________________________________________________
 364 void AliAnalysisTaskEmcalTriggerInfoQA::Terminate(Option_t *) //specify what you want to have done
 365 {
 366     // Called once at the end of the query. Done nothing here.
 367 }