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