7f3a554d3c535ba23e707b3611e9ce5730e6c838
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnMiniAnalysisTask.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 "AliRsnMiniAnalysisTask.h"
38
39 ClassImp(AliRsnMiniAnalysisTask)
40
41 //__________________________________________________________________________________________________
42 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
43    AliAnalysisTaskSE(),
44    fUseMC(kFALSE),
45    fEvNum(0),
46    fUseCentrality(kFALSE),
47    fCentralityType("QUALITY"),
48    fNMix(0),
49    fMaxDiffMult(10),
50    fMaxDiffVz(1.0),
51    fMaxDiffAngle(1E20),
52    fOutput(0x0),
53    fHistograms("AliRsnMiniOutput", 0),
54    fValues("AliRsnMiniValue", 0),
55    fHEventStat(0x0),
56    fEventCuts(0x0),
57    fTrackCuts(0),
58    fRsnEvent(),
59    fEvBuffer(0x0),
60    fTriggerAna(0x0),
61    fESDtrackCuts(0x0),
62    fMiniEvent(0x0)
63 {
64 //
65 // Dummy constructor ALWAYS needed for I/O.
66 //
67 }
68
69 //__________________________________________________________________________________________________
70 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
71    AliAnalysisTaskSE(name),
72    fUseMC(useMC),
73    fEvNum(0),
74    fUseCentrality(kFALSE),
75    fCentralityType("QUALITY"),
76    fNMix(0),
77    fMaxDiffMult(10),
78    fMaxDiffVz(1.0),
79    fMaxDiffAngle(1E20),
80    fOutput(0x0),
81    fHistograms("AliRsnMiniOutput", 0),
82    fValues("AliRsnMiniValue", 0),
83    fHEventStat(0x0),
84    fEventCuts(0x0),
85    fTrackCuts(0),
86    fRsnEvent(),
87    fEvBuffer(0x0),
88    fTriggerAna(0x0),
89    fESDtrackCuts(0x0),
90    fMiniEvent(0x0)
91 {
92 //
93 // Default constructor.
94 // Define input and output slots here (never in the dummy constructor)
95 // Input slot #0 works with a TChain - it is connected to the default input container
96 // Output slot #1 writes into a TH1 container
97 //
98
99    DefineOutput(1, TList::Class());
100 }
101
102 //__________________________________________________________________________________________________
103 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask& copy) :
104    AliAnalysisTaskSE(copy),
105    fUseMC(copy.fUseMC),
106    fEvNum(0),
107    fUseCentrality(copy.fUseCentrality),
108    fCentralityType(copy.fCentralityType),
109    fNMix(copy.fNMix),
110    fMaxDiffMult(copy.fMaxDiffMult),
111    fMaxDiffVz(copy.fMaxDiffVz),
112    fMaxDiffAngle(copy.fMaxDiffAngle),
113    fOutput(0x0),
114    fHistograms(copy.fHistograms),
115    fValues(copy.fValues),
116    fHEventStat(0x0),
117    fEventCuts(copy.fEventCuts),
118    fTrackCuts(copy.fTrackCuts),
119    fRsnEvent(),
120    fEvBuffer(0x0),
121    fTriggerAna(copy.fTriggerAna),
122    fESDtrackCuts(copy.fESDtrackCuts),
123    fMiniEvent(0x0)
124 {
125 //
126 // Copy constructor.
127 // Implemented as requested by C++ standards.
128 // Can be used in PROOF and by plugins.
129 //
130 }
131
132 //__________________________________________________________________________________________________
133 AliRsnMiniAnalysisTask& AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalysisTask& copy)
134 {
135 //
136 // Assignment operator.
137 // Implemented as requested by C++ standards.
138 // Can be used in PROOF and by plugins.
139 //
140
141    AliAnalysisTaskSE::operator=(copy);
142    
143    fUseMC = copy.fUseMC;
144    fUseCentrality = copy.fUseCentrality;
145    fCentralityType = copy.fCentralityType;
146    fNMix = copy.fNMix;
147    fMaxDiffMult = copy.fMaxDiffMult;
148    fMaxDiffVz = copy.fMaxDiffVz;
149    fMaxDiffAngle = copy.fMaxDiffAngle;
150    fHistograms = copy.fHistograms;
151    fValues = copy.fValues;
152    fEventCuts = copy.fEventCuts;
153    fTrackCuts = copy.fTrackCuts;
154    fTriggerAna = copy.fTriggerAna;
155    fESDtrackCuts = copy.fESDtrackCuts;
156    
157    return (*this);
158 }
159
160 //__________________________________________________________________________________________________
161 AliRsnMiniAnalysisTask::~AliRsnMiniAnalysisTask()
162 {
163 //
164 // Destructor. 
165 // Clean-up the output list, but not the histograms that are put inside
166 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
167 //
168
169
170    if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
171       delete fOutput;
172       delete fEvBuffer;
173    }
174 }
175
176 //__________________________________________________________________________________________________
177 Int_t AliRsnMiniAnalysisTask::AddTrackCuts(AliRsnCutSet *cuts)
178 {
179 //
180 // Add a new cut set for a new criterion for track selection.
181 // A user can add as many as he wants, and each one corresponds
182 // to one of the available bits in the AliRsnMiniParticle mask.
183 // The only check is the following: if a cut set with the same name
184 // as the argument is there, this is not added.
185 // Return value is the array position of this set.
186 //
187
188    TObject *obj = fTrackCuts.FindObject(cuts->GetName());
189    
190    if (obj) {
191       AliInfo(Form("A cut set named '%s' already exists", cuts->GetName()));
192       return fTrackCuts.IndexOf(obj);
193    } else {
194       fTrackCuts.AddLast(cuts);
195       return fTrackCuts.IndexOf(cuts);
196    }
197 }
198
199 //__________________________________________________________________________________________________
200 void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
201 {
202 //
203 // Initialization of outputs.
204 // This is called once per worker node.
205 //
206
207    // reset counter
208    fEvNum = 0;
209    
210    // message
211    AliInfo(Form("Selected event characterization: %s (%s)", (fUseCentrality ? "centrality" : "multiplicity"), fCentralityType.Data()));
212
213    // initialize trigger analysis
214    if (fTriggerAna) delete fTriggerAna;
215    fTriggerAna = new AliTriggerAnalysis;
216    
217    // initialize ESD quality cuts
218    if (fESDtrackCuts) delete fESDtrackCuts;
219    fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
220
221    // create list and set it as owner of its content (MANDATORY)
222    fOutput = new TList();
223    fOutput->SetOwner();
224    
225    // initialize event statistics counter
226    fHEventStat = new TH1F("hEventStat", "Event statistics", 4, 0.0, 4.0);
227    fHEventStat->GetXaxis()->SetBinLabel(1, "CINT1B");
228    fHEventStat->GetXaxis()->SetBinLabel(2, "V0AND");
229    fHEventStat->GetXaxis()->SetBinLabel(3, "Candle");
230    fHEventStat->GetXaxis()->SetBinLabel(4, "Accepted");
231    fOutput->Add(fHEventStat);
232    
233    // create temporary tree for filtered events
234    if (fMiniEvent) delete fMiniEvent;
235    fEvBuffer = new TTree("EventBuffer", "Temporary buffer for mini events");
236    fEvBuffer->Branch("events", "AliRsnMiniEvent", &fMiniEvent);
237    
238    // create one histogram per each stored definition (event histograms)
239    Int_t i, ndef = fHistograms.GetEntries();
240    AliRsnMiniOutput *def = 0x0;
241    for (i = 0; i < ndef; i++) {
242       def = (AliRsnMiniOutput*)fHistograms[i];
243       if (!def) continue;
244       if (!def->Init(GetName(), fOutput)) {
245          AliError(Form("Def '%s': failed initialization", def->GetName()));
246          continue;
247       }
248    }
249    
250    // post data for ALL output slots >0 here, to get at least an empty histogram
251    PostData(1, fOutput);
252 }
253
254 //__________________________________________________________________________________________________
255 void AliRsnMiniAnalysisTask::UserExec(Option_t *)
256 {
257 //
258 // Computation loop.
259 // In this case, it checks if the event is acceptable, and eventually
260 // creates the corresponding mini-event and stores it in the buffer.
261 // The real histogram filling is done at the end, in "FinishTaskOutput".
262 //
263
264    // event counter
265    fEvNum++;
266    
267    // check current event
268    Char_t check = CheckCurrentEvent();
269    if (!check) return;
270    
271    // setup PID response
272    AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager(); 
273         AliInputEventHandler *inputHandler = (AliInputEventHandler*)man->GetInputEventHandler(); 
274         fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse());
275    
276    // fill a mini-event from current
277    // and skip this event if no tracks were accepted
278    FillMiniEvent(check);
279    
280    // fill MC based histograms on mothers,
281    // which do need the original event
282    if (fUseMC) {
283       if (fRsnEvent.IsESD() && fMCEvent)
284          FillTrueMotherESD(fMiniEvent);
285       else if (fRsnEvent.IsAOD() && fRsnEvent.GetAODList())
286          FillTrueMotherAOD(fMiniEvent);
287    }
288    
289    // store event
290    fEvBuffer->Fill();
291    
292    // post data for computed stuff
293    PostData(1, fOutput);
294 }
295
296 //__________________________________________________________________________________________________
297 void AliRsnMiniAnalysisTask::FinishTaskOutput()
298 {
299 //
300 // This function is called at the end of the loop on available events,
301 // and then the buffer will be full with all the corresponding mini-events,
302 // each one containing all tracks selected by each of the available track cuts.
303 // Here a loop is done on each of these events, and both single-event and mixing are computed
304 //
305
306    // security code: reassign the buffer to the mini-event cursor
307    fEvBuffer->SetBranchAddress("events", &fMiniEvent);
308    
309    // prepare variables
310    Int_t ievt, nEvents = (Int_t)fEvBuffer->GetEntries();
311    Int_t idef, nDefs   = fHistograms.GetEntries();
312    Int_t imix, iloop, ifill;
313    AliRsnMiniOutput *def = 0x0;
314    AliRsnMiniOutput::EComputation compType;
315    
316    // loop on events, and for each one fill all outputs
317    // using the appropriate procedure depending on its type
318    // only mother-related histograms are filled in UserExec,
319    // since they require direct access to MC event
320    for (ievt = 0; ievt < nEvents; ievt++) {
321       // get next entry
322       fEvBuffer->GetEntry(ievt);
323       // fill
324       for (idef = 0; idef < nDefs; idef++) {
325          def = (AliRsnMiniOutput*)fHistograms[idef];
326          if (!def) continue;
327          compType = def->GetComputation();
328          // execute computation in the appropriate way
329          switch (compType) {
330             case AliRsnMiniOutput::kEventOnly:
331                AliDebugClass(1, Form("Event %d, def '%s': event-value histogram filling", ievt, def->GetName()));
332                ifill = 1;
333                def->FillEvent(fMiniEvent, &fValues);
334                break;
335             case AliRsnMiniOutput::kTruePair:
336                AliDebugClass(1, Form("Event %d, def '%s': true-pair histogram filling", ievt, def->GetName()));
337                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
338                break;
339             case AliRsnMiniOutput::kTrackPair:
340                AliDebugClass(1, Form("Event %d, def '%s': pair-value histogram filling", ievt, def->GetName()));
341                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
342                break;
343             case AliRsnMiniOutput::kTrackPairRotated1:
344                AliDebugClass(1, Form("Event %d, def '%s': rotated (1) background histogram filling", ievt, def->GetName()));
345                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
346                break;
347             case AliRsnMiniOutput::kTrackPairRotated2:
348                AliDebugClass(1, Form("Event %d, def '%s': rotated (2) background histogram filling", ievt, def->GetName()));
349                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
350                break;
351             default:
352                // other kinds are processed elsewhere
353                ifill = 0;
354                AliDebugClass(2, Form("Computation = %d", (Int_t)compType));
355          }
356          // message
357          AliDebugClass(1, Form("Event %6d: def = '%15s' -- fills = %5d", ievt, def->GetName(), ifill));
358       }
359    }
360    
361    // if no mixing is required, stop here and post the output
362    if (fNMix < 1) {
363       AliDebugClass(2, "Stopping here, since no mixing is required");
364       PostData(1, fOutput);
365       return;
366    }
367    
368    // initialize mixing counter
369    TString *matched = new TString[nEvents];
370    for (ievt = 0; ievt < nEvents; ievt++) {
371       matched[ievt] = "|";
372    }
373    
374    // search for good matchings
375    for (ievt = 0; ievt < nEvents; ievt++) {
376       ifill = 0;
377       fEvBuffer->GetEntry(ievt);
378       AliRsnMiniEvent evMain(*fMiniEvent);
379       for (iloop = 1; iloop < nEvents; iloop++) {
380          imix = ievt + iloop;
381          if (imix >= nEvents) imix -= nEvents;
382          if (imix == ievt) continue;
383          // text next entry
384          fEvBuffer->GetEntry(imix);
385          // skip if events are not matched
386          if (TMath::Abs(evMain.Vz()    - fMiniEvent->Vz()   ) > fMaxDiffVz   ) continue;
387          if (TMath::Abs(evMain.Mult()  - fMiniEvent->Mult() ) > fMaxDiffMult ) continue;
388          if (TMath::Abs(evMain.Angle() - fMiniEvent->Angle()) > fMaxDiffAngle) continue;
389          // check that the array of good matches for mixed does not already contain main event
390          if (matched[imix].Contains(Form("|%d|", ievt))) continue;
391          // add new mixing candidate
392          matched[ievt].Append(Form("%d|", imix));
393          ifill++;
394          if (ifill >= fNMix) break;
395       }
396       AliDebugClass(1, Form("Matches for event %5d = %s", ievt, matched[ievt].Data()));
397    }
398    
399    // perform mixing
400    TObjArray *list = 0x0;
401    TObjString *os = 0x0;
402    for (ievt = 0; ievt < nEvents; ievt++) {
403       ifill = 0;
404       fEvBuffer->GetEntry(ievt);
405       AliRsnMiniEvent evMain(*fMiniEvent);
406       list = matched[ievt].Tokenize("|");
407       TObjArrayIter next(list);
408       while ( (os = (TObjString*)next()) ) {
409          imix = os->GetString().Atoi();
410          fEvBuffer->GetEntry(imix);
411          for (idef = 0; idef < nDefs; idef++) {
412             def = (AliRsnMiniOutput*)fHistograms[idef];
413             if (!def) continue;
414             if (!def->IsTrackPairMix()) continue;
415             ifill += def->FillPair(&evMain, fMiniEvent, &fValues, kTRUE);
416             if (!def->IsSymmetric()) ifill += def->FillPair(fMiniEvent, &evMain, &fValues, kFALSE);
417          }
418       }
419       delete list;
420    }
421    
422    delete [] matched;
423       
424    /*
425    OLD   
426    ifill = 0;
427    for (iloop = 1; iloop < nEvents; iloop++) {
428       imix = ievt + iloop;
429       // restart from beginning if reached last event
430       if (imix >= nEvents) imix -= nEvents;
431       // avoid to mix an event with itself
432       if (imix == ievt) continue;
433       // skip all events already mixed enough times
434       if (fNMixed[ievt] >= fNMix) break;
435       if (fNMixed[imix] >= fNMix) continue;
436       fEvBuffer->GetEntry(imix);
437       // skip if events are not matched
438       if (TMath::Abs(evMain.Vz()    - fMiniEvent->Vz()   ) > fMaxDiffVz   ) continue;
439       if (TMath::Abs(evMain.Mult()  - fMiniEvent->Mult() ) > fMaxDiffMult ) continue;
440       if (TMath::Abs(evMain.Angle() - fMiniEvent->Angle()) > fMaxDiffAngle) continue;
441       // found a match: increment counter for both events
442       AliDebugClass(1, Form("Event %d, def '%s': event mixing (%d with %d)", ievt, def->GetName(), ievt, imix));
443       fNMixed[ievt]++;
444       fNMixed[imix]++;
445       // process mixing
446       ifill += def->FillPair(&evMain, fMiniEvent, &fValues);
447       // stop if mixed enough times
448       if (fNMixed[ievt] >= fNMix) break;
449    }
450    break;      
451    // print number of mixings done with each event
452    for (ievt = 0; ievt < nEvents; ievt++) {
453       AliDebugClass(2, Form("Event %6d: mixed %2d times", ievt, fNMixed[ievt]));
454    }
455    */
456    
457    // post computed data
458    PostData(1, fOutput);
459 }
460
461 //__________________________________________________________________________________________________
462 void AliRsnMiniAnalysisTask::Terminate(Option_t *)
463 {
464 //
465 // Draw result to screen, or perform fitting, normalizations
466 // Called once at the end of the query
467 //
468
469    fOutput = dynamic_cast<TList*>(GetOutputData(1));
470    if (!fOutput) { 
471       AliError("Could not retrieve TList fOutput"); 
472       return; 
473    }
474 }
475
476 //__________________________________________________________________________________________________
477 Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
478 {
479 //
480 // This method checks if current event is OK for analysis.
481 // In case it is, the pointers of the local AliRsnEvent data member
482 // will point to it, in order to allow cut checking, otherwise the
483 // function exits with a failure message.
484 // ---
485 // ESD events must pass the physics selection, AOD are supposed to do.
486 // ---
487 // While checking the event, a histogram is filled to count the number
488 // of CINT1B, V0AND and CANDLE events, which are needed for normalization
489 // ---
490 // Return values can be:
491 //    -- 'E' if the event is accepted and is ESD
492 //    -- 'A' if the event is accepted and is AOD
493 //    --  0  if the event is not accepted
494 //
495
496    // string to sum messages
497    TString msg("");
498    
499    // check input type
500    // exit points are provided in all cases an event is bad
501    // if this block is passed, an event can be rejected only
502    // if it does not pass the set of event cuts defined in the task
503    Char_t output = 0;
504    Bool_t isSelected;
505    if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
506       // type ESD
507       output = 'E';
508       // ESD specific check: Physics Selection
509       // --> if this is failed, the event is rejected
510       isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
511       if (!isSelected) {
512          AliDebugClass(2, "Event does not pass physics selections");
513          fRsnEvent.SetRef(0x0);
514          fRsnEvent.SetRefMC(0x0);
515          return 0;
516       }
517       // set reference to input
518       fRsnEvent.SetRef(fInputEvent);
519       // add MC if requested and available
520       if (fUseMC) {
521          if (fMCEvent) 
522             fRsnEvent.SetRefMC(fMCEvent);
523          else {
524             AliWarning("MC event requested but not available");
525             fRsnEvent.SetRefMC(0x0);
526          }
527       }
528    } else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
529       // type AOD
530       output = 'A';
531       // set reference to input
532       fRsnEvent.SetRef(fInputEvent);
533       // add MC if requested and available (it is in the same object)
534       if (fUseMC) {
535          fRsnEvent.SetRefMC(fInputEvent);
536          if (!fRsnEvent.GetAODList()) {
537             AliWarning("MC event requested but not available");
538             fRsnEvent.SetRefMC(0x0);
539          }
540       }
541    } else {
542       AliError(Form("Bad input event class: %s", fInputEvent->ClassName()));
543       // reset pointers in local AliRsnEvent object
544       fRsnEvent.SetRef(0x0);
545       fRsnEvent.SetRefMC(0x0);
546       return 0;
547    }
548    
549    // fill counter of accepted events
550    fHEventStat->Fill(0.1);
551    
552    // check if it is V0AND
553    // --> uses a cast to AliESDEvent even if the input is an AliAODEvent
554    Bool_t v0A = fTriggerAna->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0A);
555    Bool_t v0C = fTriggerAna->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0C);
556    if (v0A && v0C) {
557       msg += " -- VOAND = YES";
558       fHEventStat->Fill(1.1);
559    } else {
560       msg += " -- VOAND = NO ";
561    }
562    
563    // check candle
564    // --> requires at least one good quality track with Pt > 0.5 and |eta| <= 0.8
565    Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
566    Bool_t candle = kFALSE;
567    for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {    
568       AliVTrack   *track = (AliVTrack*)fInputEvent->GetTrack(iTrack);
569       AliESDtrack *esdt  = dynamic_cast<AliESDtrack*>(track);
570       AliAODTrack *aodt  = dynamic_cast<AliAODTrack*>(track);
571       if (track->Pt() < 0.5) continue;
572       if(TMath::Abs(track->Eta()) > 0.8) continue;
573       if (esdt) if (!fESDtrackCuts->AcceptTrack(esdt)) continue;
574       if (aodt) if (!aodt->TestFilterBit(5)) continue;
575       candle = kTRUE;
576       break;
577    }
578    if (candle) {
579       msg += " -- CANDLE = YES"; 
580       fHEventStat->Fill(2.1);
581    } else {
582       msg += " -- CANDLE = NO "; 
583    }
584    
585    // if event cuts are defined, they are checked here
586    // final decision on the event depends on this
587    isSelected = kTRUE;
588    if (fEventCuts) {
589       if (!fEventCuts->IsSelected(&fRsnEvent)) {
590          msg += " -- Local cuts = REJECTED";
591          isSelected = kFALSE;
592       } else {
593          msg += " -- Local cuts = ACCEPTED";
594          isSelected = kTRUE;
595       }
596    } else {
597       msg += " -- Local cuts = NONE";
598       isSelected = kTRUE;
599    }
600    
601    // if the above exit point is not taken, the event is accepted
602    AliDebugClass(2, Form("Stats for event %d: %s", fEvNum, msg.Data()));
603    if (isSelected) {
604       fHEventStat->Fill(3.1);
605       return output;
606    } else {
607       return 0;
608    }
609 }
610
611 //__________________________________________________________________________________________________
612 void AliRsnMiniAnalysisTask::FillMiniEvent(Char_t evType)
613 {
614 //
615 // Refresh cursor mini-event data member to fill with current event.
616 // Returns the total number of tracks selected.
617 //
618
619    // assign event-related values
620    if (fMiniEvent) delete fMiniEvent;
621    fMiniEvent = new AliRsnMiniEvent;
622    fMiniEvent->ID()    = fEvNum;
623    fMiniEvent->Vz()    = fInputEvent->GetPrimaryVertex()->GetZ();
624    fMiniEvent->Angle() = ComputeAngle();
625    fMiniEvent->Mult()  = ComputeCentrality((evType == 'E'));
626    AliDebugClass(2, Form("Event %d: type = %c -- vz = %f -- mult = %f -- angle = %f", fEvNum, evType, fMiniEvent->Vz(), fMiniEvent->Mult(), fMiniEvent->Angle()));
627    
628    // loop on daughters and assign track-related values
629    Int_t ic, ncuts = fTrackCuts.GetEntries();
630    Int_t ip, npart = fRsnEvent.GetAbsoluteSum();
631    Int_t npos = 0, nneg = 0, nneu = 0;
632    AliRsnDaughter cursor;
633    AliRsnMiniParticle miniParticle;
634    for (ip = 0; ip < npart; ip++) {
635       // point cursor to next particle
636       fRsnEvent.SetDaughter(cursor, ip);
637       // copy momentum and MC info if present
638       miniParticle.CopyDaughter(&cursor);
639       miniParticle.Index() = ip;
640       // switch on the bits corresponding to passed cuts
641       for (ic = 0; ic < ncuts; ic++) {
642          AliRsnCutSet *cuts = (AliRsnCutSet*)fTrackCuts[ic];
643          if (cuts->IsSelected(&cursor)) miniParticle.SetCutBit(ic);
644       }
645       // if a track passes at least one track cut, it is added to the pool
646       if (miniParticle.CutBits()) {
647          fMiniEvent->AddParticle(miniParticle);
648          if (miniParticle.Charge() == '+') npos++;
649          else if (miniParticle.Charge() == '-') nneg++;
650          else nneu++;
651       }
652    }
653    
654    // get number of accepted tracks
655    AliDebugClass(1, Form("Event %6d: total = %5d, accepted = %4d (pos %4d, neg %4d, neu %4d)", fEvNum, npart, (Int_t)fMiniEvent->Particles().GetEntriesFast(), npos, nneg, nneu));
656 }
657
658 //__________________________________________________________________________________________________
659 Double_t AliRsnMiniAnalysisTask::ComputeAngle()
660 {
661 //
662 // Get the plane angle
663 //
664
665    AliEventplane *plane = fInputEvent->GetEventplane();
666    if (plane) 
667       return plane->GetEventplane("Q");
668    else {
669       AliWarning("No event plane defined");
670       return 1E20;
671    }
672 }
673
674 //__________________________________________________________________________________________________
675 Double_t AliRsnMiniAnalysisTask::ComputeCentrality(Bool_t isESD)
676 {
677 //
678 // Computes event centrality/multiplicity according to the criterion defined
679 // by two elements: (1) choice between multiplicity and centrality and
680 // (2) the string defining what criterion must be used for specific computation.
681 //
682
683    if (fUseCentrality) {
684       AliCentrality *centrality = fInputEvent->GetCentrality();
685       if (!centrality) {
686          AliError("Cannot compute centrality!");
687          return -1.0;
688       }
689       return centrality->GetCentralityPercentile(fCentralityType.Data());
690    } else {
691       if (!fCentralityType.CompareTo("TRACKS"))
692          return fInputEvent->GetNumberOfTracks();
693       else if (!fCentralityType.CompareTo("QUALITY"))
694          if (isESD)
695             return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent*)fInputEvent, kTRUE);
696          else {
697             Double_t count = 0.;
698             Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
699             for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {    
700                AliVTrack   *track = (AliVTrack*)fInputEvent->GetTrack(iTrack);
701                AliAODTrack *aodt  = dynamic_cast<AliAODTrack*>(track);
702                if (!aodt) continue;
703                if (!aodt->TestFilterBit(5)) continue;
704                count++;
705             }
706             return count;
707          }
708       else if (!fCentralityType.CompareTo("TRACKLETS")) {
709          if (isESD) {
710             const AliMultiplicity *mult = ((AliESDEvent*)fInputEvent)->GetMultiplicity();
711             Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
712             for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
713             return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
714          } else {
715             AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
716             return 1E20;
717          }
718       } else {
719          AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", fCentralityType.Data()));
720          return -1.0;
721       }
722    }
723 }
724
725 //__________________________________________________________________________________________________
726 void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
727 {
728 //
729 // Fills the histograms with true mother (ESD version)
730 //
731
732    Bool_t okMatch;
733    Int_t id, ndef = fHistograms.GetEntries();
734    Int_t ip, label1, label2, npart = fMCEvent->GetNumberOfTracks();
735    static AliRsnMiniPair miniPair;
736    AliMCParticle *daughter1, *daughter2;
737    TLorentzVector p1, p2;
738    AliRsnMiniOutput *def = 0x0;
739    
740    for (id = 0; id < ndef; id++) {
741       def = (AliRsnMiniOutput*)fHistograms[id];
742       if (!def) continue;
743       if (!def->IsMother()) continue;
744       for (ip = 0; ip < npart; ip++) {
745          AliMCParticle *part = (AliMCParticle*)fMCEvent->GetTrack(ip);
746          if (TMath::Abs(part->Particle()->GetPdgCode()) != def->GetMotherPDG()) continue;
747          // check that daughters match expected species
748          if (part->Particle()->GetNDaughters() < 2) continue;
749          label1 = part->Particle()->GetDaughter(0);
750          label2 = part->Particle()->GetDaughter(1);
751          daughter1 = (AliMCParticle*)fMCEvent->GetTrack(label1);
752          daughter2 = (AliMCParticle*)fMCEvent->GetTrack(label2);
753          okMatch = kFALSE;
754          if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(1)) {
755             okMatch = kTRUE;
756             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
757             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
758          } else if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(0)) {
759             okMatch = kTRUE;
760             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
761             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
762          }
763          if (!okMatch) continue;
764          // assign momenta to computation object
765          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
766          miniPair.FillRef(def->GetMotherMass());
767          // do computations
768          def->FillMother(&miniPair, miniEvent, &fValues);
769       }
770    }
771 }
772
773 //__________________________________________________________________________________________________
774 void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
775 {
776 //
777 // Fills the histograms with true mother (AOD version)
778 //
779
780    Bool_t okMatch;
781    TClonesArray *list = fRsnEvent.GetAODList();
782    Int_t id, ndef = fHistograms.GetEntries();
783    Int_t ip, label1, label2, npart = list->GetEntries();
784    static AliRsnMiniPair miniPair;
785    AliAODMCParticle *daughter1, *daughter2;
786    TLorentzVector p1, p2;
787    AliRsnMiniOutput *def = 0x0;
788    
789    for (id = 0; id < ndef; id++) {
790       def = (AliRsnMiniOutput*)fHistograms[id];
791       if (!def) continue;
792       if (!def->IsMother()) continue;
793       for (ip = 0; ip < npart; ip++) {
794          AliAODMCParticle *part = (AliAODMCParticle*)list->At(ip);
795          if (TMath::Abs(part->GetPdgCode()) != def->GetMotherPDG()) continue;
796          // check that daughters match expected species
797          if (part->GetNDaughters() < 2) continue;
798          label1 = part->GetDaughter(0);
799          label2 = part->GetDaughter(1);
800          daughter1 = (AliAODMCParticle*)list->At(label1);
801          daughter2 = (AliAODMCParticle*)list->At(label2);
802          okMatch = kFALSE;
803          if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(1)) {
804             okMatch = kTRUE;
805             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
806             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
807          } else if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(0)) {
808             okMatch = kTRUE;
809             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
810             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
811          }
812          if (!okMatch) continue;
813          // assign momenta to computation object
814          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
815          miniPair.FillRef(def->GetMotherMass());
816          // do computations
817          def->FillMother(&miniPair, miniEvent, &fValues);
818       }
819    }
820 }