PWGLF/RESONANCES/AliRsnMiniMonitorTask.cxx

   1 //
   2 // Analysis task for 'mini' sub-package
   3 // Contains all definitions needed for running an analysis:
   4 // -- global event cut
   5 // -- a list of track cuts (any number)
   6 // -- definitions of output histograms
   7 // -- values to be computed.
   8 // Each one must be defined using the "CREATE" methods, which
   9 // add directly a new element in the task collections, and don't
  10 // need an external object to be passed to the task itself.
  11 //
  12
  13 #include <Riostream.h>
  14
  15 #include <TH1.h>
  16 #include <TList.h>
  17 #include <TTree.h>
  18
  19 #include "AliLog.h"
  20 #include "AliEventplane.h"
  21 #include "AliMultiplicity.h"
  22 #include "AliTriggerAnalysis.h"
  23 #include "AliAnalysisManager.h"
  24 #include "AliInputEventHandler.h"
  25
  26 #include "AliESDtrackCuts.h"
  27 #include "AliESDUtils.h"
  28
  29 #include "AliAODEvent.h"
  30 #include "AliAODMCParticle.h"
  31
  32 #include "AliRsnCutSet.h"
  33 #include "AliRsnMiniPair.h"
  34 #include "AliRsnMiniEvent.h"
  35 #include "AliRsnMiniParticle.h"
  36
  37 #include "AliRsnMiniMonitorTask.h"
  38
  39 ClassImp(AliRsnMiniMonitorTask)
  40
  41 //__________________________________________________________________________________________________
  42 AliRsnMiniMonitorTask::AliRsnMiniMonitorTask() :
  43    AliAnalysisTaskSE(),
  44    fUseMC(kFALSE),
  45    fEvNum(0),
  46    fUseCentrality(kFALSE),
  47    fCentralityType("QUALITY"),
  48    fOutput(0x0),
  49    fHistograms("AliRsnMiniMonitor", 0),
  50    fEventCuts(0x0),
  51    fTrackCuts(0),
  52    fRsnEvent(),
  53    fBigOutput(kFALSE)
  54 {
  55 //
  56 // Dummy constructor ALWAYS needed for I/O.
  57 //
  58 }
  59
  60 //__________________________________________________________________________________________________
  61 AliRsnMiniMonitorTask::AliRsnMiniMonitorTask(const char *name, Bool_t useMC) :
  62    AliAnalysisTaskSE(name),
  63    fUseMC(useMC),
  64    fEvNum(0),
  65    fUseCentrality(kFALSE),
  66    fCentralityType("QUALITY"),
  67    fOutput(0x0),
  68    fHistograms("AliRsnMiniMonitor", 0),
  69    fEventCuts(0x0),
  70    fTrackCuts(0),
  71    fRsnEvent(),
  72    fBigOutput(kFALSE)
  73 {
  74 //
  75 // Default constructor.
  76 // Define input and output slots here (never in the dummy constructor)
  77 // Input slot #0 works with a TChain - it is connected to the default input container
  78 // Output slot #1 writes into a TH1 container
  79 //
  80
  81    DefineOutput(1, TList::Class());
  82 }
  83
  84 //__________________________________________________________________________________________________
  85 AliRsnMiniMonitorTask::AliRsnMiniMonitorTask(const AliRsnMiniMonitorTask &copy) :
  86    AliAnalysisTaskSE(copy),
  87    fUseMC(copy.fUseMC),
  88    fEvNum(0),
  89    fUseCentrality(copy.fUseCentrality),
  90    fCentralityType(copy.fCentralityType),
  91    fOutput(0x0),
  92    fHistograms(copy.fHistograms),
  93    fEventCuts(copy.fEventCuts),
  94    fTrackCuts(copy.fTrackCuts),
  95    fRsnEvent(),
  96    fBigOutput(copy.fBigOutput)
  97 {
  98 //
  99 // Copy constructor.
 100 // Implemented as requested by C++ standards.
 101 // Can be used in PROOF and by plugins.
 102 //
 103 }
 104
 105 //__________________________________________________________________________________________________
 106 AliRsnMiniMonitorTask &AliRsnMiniMonitorTask::operator=(const AliRsnMiniMonitorTask &copy)
 107 {
 108 //
 109 // Assignment operator.
 110 // Implemented as requested by C++ standards.
 111 // Can be used in PROOF and by plugins.
 112 //
 113    AliAnalysisTaskSE::operator=(copy);
 114    if (this == &copy)
 115       return *this;
 116    fUseMC = copy.fUseMC;
 117    fUseCentrality = copy.fUseCentrality;
 118    fCentralityType = copy.fCentralityType;
 119    fHistograms = copy.fHistograms;
 120    fEventCuts = copy.fEventCuts;
 121    fTrackCuts = copy.fTrackCuts;
 122    fBigOutput = copy.fBigOutput;
 123
 124    return (*this);
 125 }
 126
 127 //__________________________________________________________________________________________________
 128 AliRsnMiniMonitorTask::~AliRsnMiniMonitorTask()
 129 {
 130 //
 131 // Destructor.
 132 // Clean-up the output list, but not the histograms that are put inside
 133 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 134 //
 135
 136
 137    if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 138       delete fOutput;
 139    }
 140 }
 141
 142 //__________________________________________________________________________________________________
 143 Int_t AliRsnMiniMonitorTask::AddTrackCuts(AliRsnCutSet *cuts)
 144 {
 145 //
 146 // Add a new cut set for a new criterion for track selection.
 147 // A user can add as many as he wants, and each one corresponds
 148 // to one of the available bits in the AliRsnMiniParticle mask.
 149 // The only check is the following: if a cut set with the same name
 150 // as the argument is there, this is not added.
 151 // Return value is the array position of this set.
 152 //
 153
 154    TObject *obj = fTrackCuts.FindObject(cuts->GetName());
 155
 156    if (obj) {
 157       AliInfo(Form("A cut set named '%s' already exists", cuts->GetName()));
 158       return fTrackCuts.IndexOf(obj);
 159    } else {
 160       fTrackCuts.AddLast(cuts);
 161       return fTrackCuts.IndexOf(cuts);
 162    }
 163 }
 164
 165 //__________________________________________________________________________________________________
 166 void AliRsnMiniMonitorTask::UserCreateOutputObjects()
 167 {
 168 //
 169 // Initialization of outputs.
 170 // This is called once per worker node.
 171 //
 172
 173    // reset counter
 174    fEvNum = 0;
 175
 176    // message
 177    AliInfo(Form("Selected event characterization: %s (%s)", (fUseCentrality ? "centrality" : "multiplicity"), fCentralityType.Data()));
 178
 179    // create list and set it as owner of its content (MANDATORY)
 180    if (fBigOutput) OpenFile(1);
 181    fOutput = new TList();
 182    fOutput->SetOwner();
 183
 184    // create one histogram per each stored definition (event histograms)
 185    Int_t i, ndef = fHistograms.GetEntries();
 186    AliRsnMiniMonitor *def = 0x0;
 187    for (i = 0; i < ndef; i++) {
 188       def = (AliRsnMiniMonitor *)fHistograms[i];
 189       if (!def) continue;
 190       if (!def->Init(GetName(), fOutput)) {
 191          AliError(Form("Def '%s': failed initialization", def->GetName()));
 192          continue;
 193       }
 194    }
 195
 196    // post data for ALL output slots >0 here, to get at least an empty histogram
 197    PostData(1, fOutput);
 198 }
 199
 200 //__________________________________________________________________________________________________
 201 void AliRsnMiniMonitorTask::UserExec(Option_t *)
 202 {
 203 //
 204 // Computation loop.
 205 // In this case, it checks if the event is acceptable, and eventually
 206 // creates the corresponding mini-event and stores it in the buffer.
 207 // The real histogram filling is done at the end, in "FinishTaskOutput".
 208 //
 209
 210    // event counter
 211    fEvNum++;
 212
 213    // check current event
 214    Char_t check = CheckCurrentEvent();
 215    if (!check) return;
 216
 217    // setup PID response
 218    AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 219    AliInputEventHandler *inputHandler = (AliInputEventHandler *)man->GetInputEventHandler();
 220    fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse());
 221
 222    // loop on monitors and fill them
 223    Int_t it, icut, nTracks = fRsnEvent.GetAbsoluteSum(), nCuts = fTrackCuts.GetEntriesFast();
 224    AliRsnDaughter cursor;
 225    AliRsnCutSet *cut = 0x0;
 226    AliRsnMiniMonitor *mon = 0x0;
 227    TObjArrayIter next(&fHistograms);
 228    for (it = 0; it < nTracks; it++) {
 229       fRsnEvent.SetDaughter(cursor, it, fUseMC);
 230       next.Reset();
 231       while ( (mon = (AliRsnMiniMonitor *)next()) ) {
 232          icut = mon->GetCutID();
 233          if (icut >= 0 && icut < nCuts) {
 234             cut = (AliRsnCutSet *)fTrackCuts[icut];
 235             if (!cut) {
 236                AliError("Cut not found");
 237                continue;
 238             }
 239             if (!cut->IsSelected(&cursor)) continue;
 240          }
 241          mon->Fill(&cursor, &fRsnEvent);
 242       }
 243    }
 244
 245    // post data for computed stuff
 246    PostData(1, fOutput);
 247 }
 248
 249 //__________________________________________________________________________________________________
 250 void AliRsnMiniMonitorTask::Terminate(Option_t *)
 251 {
 252 //
 253 // Draw result to screen, or perform fitting, normalizations
 254 // Called once at the end of the query
 255 //
 256
 257    fOutput = dynamic_cast<TList *>(GetOutputData(1));
 258    if (!fOutput) {
 259       AliError("Could not retrieve TList fOutput");
 260       return;
 261    }
 262 }
 263
 264 //__________________________________________________________________________________________________
 265 Char_t AliRsnMiniMonitorTask::CheckCurrentEvent()
 266 {
 267 //
 268 // This method checks if current event is OK for analysis.
 269 // In case it is, the pointers of the local AliRsnEvent data member
 270 // will point to it, in order to allow cut checking, otherwise the
 271 // function exits with a failure message.
 272 // ---
 273 // ESD events must pass the physics selection, AOD are supposed to do.
 274 // ---
 275 // While checking the event, a histogram is filled to count the number
 276 // of CINT1B, V0AND and CANDLE events, which are needed for normalization
 277 // ---
 278 // Return values can be:
 279 //    -- 'E' if the event is accepted and is ESD
 280 //    -- 'A' if the event is accepted and is AOD
 281 //    --  0  if the event is not accepted
 282 //
 283
 284    // string to sum messages
 285    TString msg("");
 286
 287    // check input type
 288    // exit points are provided in all cases an event is bad
 289    // if this block is passed, an event can be rejected only
 290    // if it does not pass the set of event cuts defined in the task
 291    Char_t output = 0;
 292    Bool_t isSelected;
 293    if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
 294       // type ESD
 295       output = 'E';
 296       // ESD specific check: Physics Selection
 297       // --> if this is failed, the event is rejected
 298       isSelected = (((AliInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 299       if (!isSelected) {
 300          AliDebugClass(2, "Event does not pass physics selections");
 301          fRsnEvent.SetRef(0x0);
 302          fRsnEvent.SetRefMC(0x0);
 303          return 0;
 304       }
 305       // set reference to input
 306       fRsnEvent.SetRef(fInputEvent);
 307       // add MC if requested and available
 308       if (fUseMC) {
 309          if (fMCEvent)
 310             fRsnEvent.SetRefMC(fMCEvent);
 311          else {
 312             AliWarning("MC event requested but not available");
 313             fRsnEvent.SetRefMC(0x0);
 314          }
 315       }
 316    } else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
 317       // type AOD
 318       output = 'A';
 319       // set reference to input
 320       fRsnEvent.SetRef(fInputEvent);
 321       // add MC if requested and available (it is in the same object)
 322       if (fUseMC) {
 323          fRsnEvent.SetRefMC(fInputEvent);
 324          if (!fRsnEvent.GetAODList()) {
 325             AliWarning("MC event requested but not available");
 326             fRsnEvent.SetRefMC(0x0);
 327          }
 328       }
 329    } else {
 330       AliError(Form("Bad input event class: %s", fInputEvent->ClassName()));
 331       // reset pointers in local AliRsnEvent object
 332       fRsnEvent.SetRef(0x0);
 333       fRsnEvent.SetRefMC(0x0);
 334       return 0;
 335    }
 336
 337    // if event cuts are defined, they are checked here
 338    // final decision on the event depends on this
 339    isSelected = kTRUE;
 340    if (fEventCuts) {
 341       if (!fEventCuts->IsSelected(&fRsnEvent)) {
 342          msg += " -- Local cuts = REJECTED";
 343          isSelected = kFALSE;
 344       } else {
 345          msg += " -- Local cuts = ACCEPTED";
 346          isSelected = kTRUE;
 347       }
 348    } else {
 349       msg += " -- Local cuts = NONE";
 350       isSelected = kTRUE;
 351    }
 352
 353    // if the above exit point is not taken, the event is accepted
 354    AliDebugClass(2, Form("Stats for event %d: %s", fEvNum, msg.Data()));
 355    if (isSelected) {
 356       return output;
 357    } else {
 358       return 0;
 359    }
 360 }
 361
 362 //__________________________________________________________________________________________________
 363 Double_t AliRsnMiniMonitorTask::ComputeCentrality(Bool_t isESD)
 364 {
 365 //
 366 // Computes event centrality/multiplicity according to the criterion defined
 367 // by two elements: (1) choice between multiplicity and centrality and
 368 // (2) the string defining what criterion must be used for specific computation.
 369 //
 370
 371    if (fUseCentrality) {
 372       AliCentrality *centrality = fInputEvent->GetCentrality();
 373       if (!centrality) {
 374          AliError("Cannot compute centrality!");
 375          return -1.0;
 376       }
 377       return centrality->GetCentralityPercentile(fCentralityType.Data());
 378    } else {
 379       if (!fCentralityType.CompareTo("TRACKS"))
 380          return fInputEvent->GetNumberOfTracks();
 381       else if (!fCentralityType.CompareTo("QUALITY"))
 382          if (isESD)
 383             return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
 384          else {
 385             Double_t count = 0.;
 386             Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
 387             for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
 388                AliVTrack   *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
 389                AliAODTrack *aodt  = dynamic_cast<AliAODTrack *>(track);
 390                if (!aodt) continue;
 391                if (!aodt->TestFilterBit(5)) continue;
 392                count++;
 393             }
 394             return count;
 395          }
 396       else if (!fCentralityType.CompareTo("TRACKLETS")) {
 397          if (isESD) {
 398             const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
 399             Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
 400             for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
 401             return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
 402          } else {
 403             AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
 404             return 1E20;
 405          }
 406       } else {
 407          AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", fCentralityType.Data()));
 408          return -1.0;
 409       }
 410    }
 411 }