346595e7787cf895ab188cbebc7776c2aa8b7e1d
[u/mrichter/AliRoot.git] / PWGCF / Correlations / Base / AliTwoPlusOneContainer.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //
17 //
18 // encapsulates several AliUEHist objects for a full UE analysis plus additional control histograms
19 // the method FillCorrelations is the only method which fills the container
20 // associated particles are filled once for each near side and once for each away side
21 // the pT 2 for the near side particles is set to the highest pT of all accepted away side pT 2
22 // as consequence the pT 2 axis can only be cut for the minimum
23 // if the maximum of the axis should be changed the full analysis needs to rerun
24 //
25 //
26 // Author: Markus Zimmermann
27 // 
28
29
30 #include "AliTwoPlusOneContainer.h"
31
32 #include "TList.h"
33
34 #include "AliCFContainer.h"
35 #include "AliVParticle.h"
36
37 #include "TH1F.h"
38 #include "TH2F.h"
39 #include "TMath.h"
40
41 ClassImp(AliTwoPlusOneContainer)
42
43 AliTwoPlusOneContainer::AliTwoPlusOneContainer(const char* name, const char* binning, Double_t alpha) : 
44   TNamed(name, name),
45   fTwoPlusOne(0),
46   fAsymmetry(0),
47   fAsymmetryMixed(0),
48   fTriggerPt(0),
49   fTriggerPt1Min(0),
50   fTriggerPt1Max(0),
51   fTriggerPt2Min(0),
52   fTriggerPt2Max(0),
53   fPtAssocMin(0),
54   fPtAssocMax(0),
55   fAlpha(alpha),
56   fUseLeadingPt(1),
57   fMergeCount(0)
58 {
59   // Constructor
60   //
61   TString defaultBinningStr;
62   defaultBinningStr = "eta: -1.0, 1.0\n"//needs to be defined for AliUEHist
63     "p_t_assoc: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0\n"
64     "p_t_leading: 6.0, 8.0, 10.0, 12.0, 14.0\n"
65     "p_t_leading_course: 4.0, 6.0, 8.0, 10.0\n"
66     "p_t_eff: 0.5, 8.0\n"//needs to be defined for AliUEHist
67     "vertex_eff: -10, 10\n"//needs to be defined for AliUEHist
68     "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n"
69     "vertex: -7, -5, -3, -1, 1, 3, 5, 7\n"
70     "delta_phi: -1.570796, -1.396263, -1.221730, -1.047197, -0.872664, -0.698131, -0.523598, -0.349065, -0.174532, 0.0, 0.174532, 0.349065, 0.523598, 0.698131, 0.872664, 1.047197, 1.221730, 1.396263, 1.570796, 1.745329, 1.919862, 2.094395, 2.268928, 2.443460, 2.617993, 2.792526, 2.967059, 3.141592, 3.316125, 3.490658, 3.665191, 3.839724, 4.014257, 4.188790, 4.363323, 4.537856, 4.712388\n" // this binning starts at -pi/2 and it ends at +3pi/2, it has 36 bins
71     "delta_eta: -1.8, -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8\n";
72
73   // combine customBinning with defaultBinningStr -> use customBinning where available and otherwise defaultBinningStr
74   TString binningStr = AliUEHist::CombineBinning(defaultBinningStr, TString(binning));
75
76   fTwoPlusOne = new AliUEHist("TwoPlusOne", binningStr);
77
78   //set minimum and maximum trigger pt values
79   fTriggerPt1Min = fTwoPlusOne->GetTrackHist(AliUEHist::kToward)->GetGrid((AliUEHist::CFStep) AliTwoPlusOneContainer::kSameNS)->GetAxis(2)->GetXmin();
80   fTriggerPt1Max = fTwoPlusOne->GetTrackHist(AliUEHist::kToward)->GetGrid((AliUEHist::CFStep) AliTwoPlusOneContainer::kSameNS)->GetAxis(2)->GetXmax();
81   fTriggerPt2Min = fTwoPlusOne->GetTrackHist(AliUEHist::kToward)->GetGrid((AliUEHist::CFStep) AliTwoPlusOneContainer::kSameNS)->GetAxis(6)->GetXmin();
82   fTriggerPt2Max = fTwoPlusOne->GetTrackHist(AliUEHist::kToward)->GetGrid((AliUEHist::CFStep) AliTwoPlusOneContainer::kSameNS)->GetAxis(6)->GetXmax();
83   fPtAssocMin = fTwoPlusOne->GetTrackHist(AliUEHist::kToward)->GetGrid((AliUEHist::CFStep) AliTwoPlusOneContainer::kSameNS)->GetAxis(1)->GetXmin();
84   fPtAssocMax = fTwoPlusOne->GetTrackHist(AliUEHist::kToward)->GetGrid((AliUEHist::CFStep) AliTwoPlusOneContainer::kSameNS)->GetAxis(1)->GetXmax();
85   
86   fAsymmetry  = new TH1F("fAsymmetry", ";A;dN", 50, 0, 1);
87   fAsymmetryMixed  = new TH1F("fAsymmetryMixed", ";A;dN", 50, 0, 1);
88   Int_t pt1_bins = (fTriggerPt1Max - fTriggerPt1Min)/0.1;
89   Int_t pt2_bins = (fTriggerPt2Max - fTriggerPt2Min)/0.1;
90   fTriggerPt = new TH2F("fTriggerPt", ";p_{T,1};p_{T,2}", pt1_bins, fTriggerPt1Min, fTriggerPt1Max, pt2_bins, fTriggerPt2Min, fTriggerPt2Max);
91   
92   TH1::AddDirectory();
93 }
94
95 //_____________________________________________________________________________
96 AliTwoPlusOneContainer::AliTwoPlusOneContainer(const AliTwoPlusOneContainer &c) :
97   TNamed(fName, fTitle),
98   fTwoPlusOne(0),
99   fAsymmetry(0),
100   fAsymmetryMixed(0),
101   fTriggerPt(0),
102   fTriggerPt1Min(0),
103   fTriggerPt1Max(0),
104   fTriggerPt2Min(0),
105   fTriggerPt2Max(0),
106   fPtAssocMin(0),
107   fPtAssocMax(0),
108   fAlpha(0.2),
109   fMergeCount(0)
110 {
111   //
112   // AliTwoPlusOneContainer copy constructor
113   //
114
115   ((AliTwoPlusOneContainer &) c).Copy(*this);
116   fTriggerPt1Min = ((AliTwoPlusOneContainer &) c).getTriggerPt1Min();
117   fTriggerPt1Max = ((AliTwoPlusOneContainer &) c).getTriggerPt1Max();
118   fTriggerPt2Min = ((AliTwoPlusOneContainer &) c).getTriggerPt2Min();
119   fTriggerPt2Max = ((AliTwoPlusOneContainer &) c).getTriggerPt2Max();
120   fPtAssocMin = ((AliTwoPlusOneContainer &) c).getPtAssocMin();
121   fPtAssocMax = ((AliTwoPlusOneContainer &) c).getPtAssocMax();
122   fAlpha = ((AliTwoPlusOneContainer &) c).fAlpha;
123   fMergeCount = ((AliTwoPlusOneContainer &) c).fMergeCount;
124 }
125
126 //____________________________________________________________________
127 AliTwoPlusOneContainer::~AliTwoPlusOneContainer()
128 {
129   // Destructor
130   
131   DeleteContainers();
132 }
133
134 //____________________________________________________________________
135 void AliTwoPlusOneContainer::DeleteContainers()
136 {
137   if (fTwoPlusOne)
138   {
139     delete fTwoPlusOne;
140     fTwoPlusOne = 0;
141   }
142   
143   if(fAsymmetry)
144   {
145     delete fAsymmetry;
146     fAsymmetry = 0;
147   }
148
149   if(fAsymmetryMixed)
150   {
151     delete fAsymmetryMixed;
152     fAsymmetryMixed = 0;
153   }
154
155   if(fTriggerPt)
156   {
157     delete fTriggerPt;
158     fTriggerPt = 0;
159   }
160 }
161
162
163 //____________________________________________________________________
164 void AliTwoPlusOneContainer::FillCorrelations(Double_t centrality, Float_t zVtx, AliTwoPlusOneContainer::PlotKind step, TObjArray* triggerNear, TObjArray* triggerAway, TObjArray* assocNear, TObjArray* assocAway, Double_t weight, Bool_t is1plus1, Bool_t isBackgroundSame)
165 {
166   //Fill Correlations fills the UEHist fTwoPlusOne with the 2+1 correlation
167   //the input variables centrality and zVtx are the centrality and the z vertex of the event
168   //the 4 lists triggerNear, triggerAway, assocNear and assocAway are four lists of particles. For the same event analysis all for lists are identical. For the mixed event analysis assocNear and assocAway are from different events as the trigger lists and for the mixed combinatorics triggerNear and assocNear are from one event and triggerAway and assocAway are from another event.
169
170   AliCFContainer* track_hist = fTwoPlusOne->GetTrackHist(AliUEHist::kToward);
171   AliCFContainer* event_hist = fTwoPlusOne->GetEventHist();
172   AliUEHist::CFStep stepUEHist = static_cast<AliUEHist::CFStep>(step);
173
174   //in case of the computation of the background in the same event there are two possible positions: delta phi = +/- pi/2
175   //both positions are used so the results could only be weighted with 0.5*weight
176   if(isBackgroundSame)
177     weight *= 0.5;
178
179   for (Int_t i=0; i<triggerNear->GetEntriesFast(); i++){
180     AliVParticle* part = (AliVParticle*) triggerNear->UncheckedAt(i);
181     
182     Double_t part_pt = part->Pt();
183     if(part_pt<fTriggerPt1Min || part_pt>fTriggerPt1Max)
184       continue;
185
186     Double_t part_eta = part->Eta();
187     Double_t part_phi = part->Phi();
188
189     //search for second trigger particle
190     //entries is the number of particles in the event, so there can only be less particles within the trigger 2 conditions
191     Int_t ind_found = 0;
192     Int_t ind_max_found_pt = -1;
193     Int_t triggerAway_entries = triggerAway->GetEntriesFast();
194     AliVParticle* found_particle[triggerAway_entries];
195
196     Bool_t do_not_use_T1 = false;
197
198     //have to fake the away side triggers for the 1+1 analysis
199     if(is1plus1){
200       found_particle[ind_found] = part;//in 1plus1 use first trigger particle also as pseudo second trigger particle
201       ind_max_found_pt = ind_found;
202       ind_found = 1;
203     }else{
204       //normal 2+1 analysis
205       for (Int_t j=0; j<triggerAway_entries; j++){
206         AliVParticle* part2 = (AliVParticle*) triggerAway->UncheckedAt(j);
207
208         Double_t part2_pt = part2->Pt();
209         //check if pT of trigger 2 has enough energy to be a trigger
210         //maximum energy is checked later after checking this particle may have to much energy for trigger 1 or 2
211         if(part2_pt<fTriggerPt2Min)
212           continue;
213
214         // don't use the same particle (is in any case impossible because the Delta phi angle will be 0)
215         if(part==part2){
216           continue;
217         }
218         
219         Double_t dphi_triggers = part_phi-part2->Phi();
220         if(dphi_triggers>1.5*TMath::Pi()) dphi_triggers -= TMath::TwoPi();
221         else if(dphi_triggers<-0.5*TMath::Pi()) dphi_triggers += TMath::TwoPi();
222       
223         if(fUseLeadingPt && part2_pt>part_pt){
224           if(TMath::Abs(dphi_triggers)<fAlpha){
225             do_not_use_T1 = true;
226             break;
227           }
228         }
229
230         //if 2+1 analysis check if trigger particles have a delta phi = pi +/- alpha
231         if(!isBackgroundSame)
232           dphi_triggers -= TMath::Pi();
233         else if(isBackgroundSame){
234           //shift defined area of delta phi
235           if(dphi_triggers>TMath::Pi()) dphi_triggers -= TMath::TwoPi();
236         
237           //look at delta phi = +/- pi/2
238           if(dphi_triggers<0)
239             dphi_triggers += 0.5*TMath::Pi();
240           else if(dphi_triggers>0)
241             dphi_triggers -= 0.5*TMath::Pi();
242         }
243         if(TMath::Abs(dphi_triggers)>fAlpha)
244           continue;
245
246         //check if pT of trigger 2 is too high
247         if(part2_pt>fTriggerPt2Max || part2_pt>part_pt){
248           //also pt of trigger 2 needs to be smaller than the pt of trigger 1 (to have an ordering if both pt are close to each other)
249           if(fUseLeadingPt){
250             do_not_use_T1 = true;
251             break;
252           }else
253             continue;
254         }
255
256         found_particle[ind_found] = part2;
257         if(ind_max_found_pt==-1 || part2_pt>found_particle[ind_max_found_pt]->Pt()) ind_max_found_pt = ind_found;
258         ind_found++;
259       }//end loop to search for the second trigger particle
260     }
261
262     //if there is a particle with higher energy than T1 closer than alpha to T1, do not use this T1
263     //if there is a particle with higher energy than T1 or max(T2) within Delta phi = pi +/- alpha to T1, do not use this T1
264     if(do_not_use_T1)
265       continue;
266
267     //if no second trigger particle was found continue to search for the next first trigger particle
268     if(ind_found==0)
269       continue;
270
271     //the energy of the second trigger particle is set for the near side to the maximum energy of all trigger 2 particles on the away side
272     // this leads to the fact that the number of accepted trigger combinations can be artificial smaller than the real number if there is a cut on the pT 2 energy from the top; cutting away the smallest energy of pT 2 is still save; this is the reason why it is not allowed to use a cut on the top pt of trigger particle 2
273     //fill trigger particles
274     if(ind_found>0){
275         Double_t vars[4];
276         vars[0] = part_pt;
277         vars[1] = centrality;
278         vars[2] = zVtx;
279         vars[3] = found_particle[ind_max_found_pt]->Pt();
280         if(is1plus1)
281           vars[3] = (fTriggerPt2Max+fTriggerPt2Min)/2;
282
283         event_hist->Fill(vars, stepUEHist, weight);//near side
284
285         if(!is1plus1)
286           for(Int_t k=0; k< ind_found; k++){
287             vars[3] = found_particle[k]->Pt();
288             event_hist->Fill(vars, stepUEHist+1, weight);//away side
289           }
290         
291         //fill fTriggerPt only once, choosed kSameNS
292         if(step==AliTwoPlusOneContainer::kSameNS)
293           for(Int_t k=0; k< ind_found; k++)
294             fTriggerPt->Fill(part_pt, found_particle[k]->Pt());
295
296         //fill asymmetry only for kSameNS and kMixedNS
297         if(step==AliTwoPlusOneContainer::kSameNS||step==AliTwoPlusOneContainer::kMixedNS){
298           for(Int_t k=0; k< ind_found; k++){
299             Float_t asymmetry = (part_pt-found_particle[k]->Pt())/(part_pt+found_particle[k]->Pt());
300             if(step==AliTwoPlusOneContainer::kSameNS){
301               fAsymmetry->Fill(asymmetry);
302             }else{
303               fAsymmetryMixed->Fill(asymmetry);
304             }
305           }
306         }
307           
308     }
309  
310     //add correlated particles on the near side
311     for (Int_t k=0; k<assocNear->GetEntriesFast(); k++){
312       AliVParticle* part3 = (AliVParticle*) assocNear->UncheckedAt(k);
313
314       Double_t part3_pt = part3->Pt();
315       if(part3_pt<fPtAssocMin || part3_pt>fPtAssocMax)
316         continue;
317
318       //do not add the trigger 1 particle
319       if(part==part3)
320         continue;
321
322       //do not add the trigger 2 particle
323       if(found_particle[ind_max_found_pt]==part3)
324         continue;
325
326       Double_t dphi_near = part_phi-part3->Phi(); 
327       if(dphi_near>1.5*TMath::Pi()) dphi_near -= TMath::TwoPi();
328       else if(dphi_near<-0.5*TMath::Pi()) dphi_near += TMath::TwoPi();
329
330       Double_t deta_near = part_eta-part3->Eta();
331       
332       Double_t vars[7];
333       vars[0] = deta_near;
334       vars[1] = part3_pt;
335       vars[2] = part_pt;
336       vars[3] = centrality;
337       vars[4] = dphi_near;
338       vars[5] = zVtx;
339       vars[6] = found_particle[ind_max_found_pt]->Pt();
340       if(is1plus1)
341         vars[6] = (fTriggerPt2Max+fTriggerPt2Min)/2;
342
343       track_hist->Fill(vars, stepUEHist, weight);
344     }
345
346     //search only for the distribution of the 2nd trigger particle
347     if(is1plus1)
348       continue;
349
350     //add correlated particles on the away side
351     for (Int_t k=0; k<assocAway->GetEntriesFast(); k++){
352       AliVParticle* part3 = (AliVParticle*) assocAway->UncheckedAt(k);
353         
354       Double_t part3_pt = part3->Pt();
355       if(part3_pt<fPtAssocMin || part3_pt>fPtAssocMax)
356         continue;
357
358       //do not add the trigger 1 particle
359       if(part==part3)
360         continue;
361
362       for(int l=0; l<ind_found; l++){
363         //do not add the trigger 2 particle
364         if(found_particle[l]==part3)
365           continue;
366
367         Double_t dphi_away = found_particle[l]->Phi()-part3->Phi();
368         if(dphi_away>1.5*TMath::Pi()) dphi_away -= TMath::TwoPi();
369         else if(dphi_away<-0.5*TMath::Pi()) dphi_away += TMath::TwoPi();
370         
371         Double_t deta_away = found_particle[l]->Eta()-part3->Eta();
372       
373         Double_t vars[7];
374         vars[0] = deta_away;
375         vars[1] = part3_pt;
376         vars[2] = part_pt;
377         vars[3] = centrality;
378         vars[4] = dphi_away;
379         vars[5] = zVtx;
380         vars[6] = found_particle[l]->Pt();
381         
382         track_hist->Fill(vars, stepUEHist+1, weight);//step +1 is the AS to the NS plot of step
383       }
384     }
385   }//end loop to search for the first trigger particle
386 }
387
388
389
390 //____________________________________________________________________
391 AliTwoPlusOneContainer &AliTwoPlusOneContainer::operator=(const AliTwoPlusOneContainer &c)
392 {
393   // assigment operator
394
395   DeleteContainers();
396
397   if (this != &c)
398     ((AliTwoPlusOneContainer &) c).Copy(*this);
399
400   return *this;
401 }
402
403 //____________________________________________________________________
404 void AliTwoPlusOneContainer::Copy(TObject& c) const
405 {
406   // copy function
407
408   AliTwoPlusOneContainer& target = (AliTwoPlusOneContainer &) c;
409
410   if (fTwoPlusOne)
411     target.fTwoPlusOne = dynamic_cast<AliUEHist*> (fTwoPlusOne->Clone());
412   
413   if (fAsymmetry)
414     target.fAsymmetry = dynamic_cast<TH1F*> (fAsymmetry->Clone());
415
416   if (fAsymmetryMixed)
417     target.fAsymmetryMixed = dynamic_cast<TH1F*> (fAsymmetryMixed->Clone());
418
419   if (fTriggerPt)
420     target.fTriggerPt = dynamic_cast<TH2F*> (fTriggerPt->Clone());
421   
422 }
423
424 //____________________________________________________________________
425 Long64_t AliTwoPlusOneContainer::Merge(TCollection* list)
426 {
427   // Merge a list of AliTwoPlusOneContainer objects with this. 
428   // Returns the number of merged objects (including this).
429
430   if (!list)
431     return 0;
432
433   if (list->IsEmpty())
434     return 1;
435
436   TIterator* iter = list->MakeIterator();
437   TObject* obj;
438   
439   // collections of objects
440   const Int_t kMaxLists = 4;
441   TList* lists[kMaxLists];
442
443   for (Int_t i=0; i<kMaxLists; i++)
444     lists[i] = new TList;
445
446   Int_t count = 0;
447   while ((obj = iter->Next())) {
448     
449     AliTwoPlusOneContainer* entry = dynamic_cast<AliTwoPlusOneContainer*> (obj);
450     if (entry == 0) 
451       continue;
452
453     lists[0]->Add(entry->fTwoPlusOne);
454     lists[1]->Add(entry->fAsymmetry);
455     lists[2]->Add(entry->fAsymmetryMixed);
456     lists[3]->Add(entry->fTriggerPt);
457
458     fMergeCount += entry->fMergeCount;
459     count++;
460   }
461   
462   fTwoPlusOne->Merge(lists[0]);
463   fAsymmetry->Merge(lists[1]);
464   fAsymmetryMixed->Merge(lists[2]);
465   fTriggerPt->Merge(lists[3]);
466
467   for (Int_t i=0; i<kMaxLists; i++)
468   delete lists[i];
469
470   //  delete lists;
471   return count+1;
472
473 }