switch off shriking of tclones array to avoid large virtual memory
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronMixingHandler.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2009, 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 //                Dielectron MixingHandler                                  //
18 //                                                                       //
19 //                                                                       //
20 /*
21 Detailed description
22
23
24 */
25 //                                                                       //
26 ///////////////////////////////////////////////////////////////////////////
27
28 #include <TVectorD.h>
29 #include <TH1.h>
30 #include <TAxis.h>
31
32 #include <AliLog.h>
33 #include <AliVTrack.h>
34
35 #include "AliDielectron.h"
36 #include "AliDielectronHelper.h"
37 #include "AliDielectronHistos.h"
38 #include "AliDielectronEvent.h"
39
40 #include "AliDielectronMixingHandler.h"
41
42 ClassImp(AliDielectronMixingHandler)
43
44 AliDielectronMixingHandler::AliDielectronMixingHandler() :
45   TNamed(),
46   fDepth(10),
47   fArrPools("TClonesArray"),
48   fAxes(kMaxCuts),
49   fMixType(kOSonly),
50   fMixIncomplete(kTRUE),
51   fMoveToSameVertex(kFALSE),
52   fSkipFirstEvt(kFALSE),
53   fPID(0x0)
54 {
55   //
56   // Default Constructor
57   //
58   for (Int_t i=0; i<kMaxCuts; ++i){
59     fEventCuts[i]=0;
60   }
61   fAxes.SetOwner(kTRUE);
62 }
63
64 //______________________________________________
65 AliDielectronMixingHandler::AliDielectronMixingHandler(const char* name, const char* title) :
66   TNamed(name, title),
67   fDepth(10),
68   fArrPools("TClonesArray"),
69   fAxes(kMaxCuts),
70   fMixType(kOSonly),
71   fMixIncomplete(kTRUE),
72   fMoveToSameVertex(kFALSE),
73   fSkipFirstEvt(kFALSE),
74   fPID(0x0)
75 {
76   //
77   // Named Constructor
78   //
79   for (Int_t i=0; i<kMaxCuts; ++i){
80     fEventCuts[i]=0;
81   }
82   fAxes.SetOwner(kTRUE);
83 }
84
85 //______________________________________________
86 AliDielectronMixingHandler::~AliDielectronMixingHandler()
87 {
88   //
89   // Default Destructor
90   //
91   fAxes.Delete();
92   delete fPID;
93 }
94
95 //________________________________________________________________
96 void AliDielectronMixingHandler::AddVariable(AliDielectronVarManager::ValueTypes type,
97                                              Int_t nbins, Double_t min, Double_t max, Bool_t log)
98 {
99   //
100   // Add a variable to the mixing handler
101   //
102
103   // limit number of variables to kMaxCuts
104   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
105   
106   TVectorD *binLimits=0x0;
107   if (!log) binLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
108   else binLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
109   if (!binLimits) return;
110
111   Int_t size=fAxes.GetEntriesFast();
112   fEventCuts[size]=(UShort_t)type;
113   fAxes.Add(binLimits);
114 }
115
116 //________________________________________________________________
117 void AliDielectronMixingHandler::AddVariable(AliDielectronVarManager::ValueTypes type,
118                                              const char* binLimitStr)
119 {
120   //
121   // Add a variable to the mixing handler with arbitrary binning
122   //
123
124   // limit number of variables to kMaxCuts
125   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
126   
127   TVectorD *binLimits=AliDielectronHelper::MakeArbitraryBinning(binLimitStr);
128   if (!binLimits) return;
129   
130   Int_t size=fAxes.GetEntriesFast();
131   fEventCuts[size]=(UShort_t)type;
132   fAxes.Add(binLimits);
133 }
134
135 //______________________________________________
136 void AliDielectronMixingHandler::Fill(const AliVEvent *ev, AliDielectron *diele)
137 {
138   //
139   // fill event buffers and perform mixing if the pool depth is reached
140   //
141
142   //check if there are tracks available
143   if (diele->GetTrackArray(0)->GetEntriesFast()==0 && diele->GetTrackArray(1)->GetEntriesFast()==0) return;
144
145   TString dim;
146   Int_t bin=FindBin(AliDielectronVarManager::GetData(),&dim);
147
148   //add mixing bin to event data
149   AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
150
151   if (bin<0){
152     AliDebug(5,Form("Bin outside range: %s",dim.Data()));
153     return;
154   }
155
156   // get mixing pool, create it if it does not yet exist.
157   TClonesArray *poolp=static_cast<TClonesArray*>(fArrPools.At(bin));
158   if (!poolp){
159     AliDebug(10,Form("New pool at %d (%s)\n",bin,dim.Data()));
160     poolp=new(fArrPools[bin]) TClonesArray("AliDielectronEvent",1);
161   }
162   TClonesArray &pool=*poolp;
163
164   // clear the current pool if its size was reached by last event
165   // clear before fill new event into it
166   // NOTE: clear not directly after DoMixing, because you may want to use the ME
167   // in the internal train by other configs/tasks
168   // reset the event pool size to 1 (this keeps the physical memory consumption low)
169   if(pool.GetEntriesFast()==fDepth)  {
170     pool.Clear("C");
171     //    pool.ExpandCreate(1);
172   }
173
174   AliDebug(10,Form("new event at %d: %d",bin,pool.GetEntriesFast()));
175   AliDielectronEvent *event=new(pool[pool.GetEntriesFast()]) AliDielectronEvent();
176   if(ev->IsA() == AliAODEvent::Class()) event->SetAOD(TMath::Max(diele->GetTrackArray(0)->GetEntriesFast(),diele->GetTrackArray(1)->GetEntriesFast()));
177   else event->SetESD();
178
179   event->SetProcessID(fPID);
180   event->SetTracks(*diele->GetTrackArray(0), *diele->GetTrackArray(1), *diele->GetPairArray(1));
181   event->SetEventData(AliDielectronVarManager::GetData());
182
183   // check if pool depth is reached.
184   if (pool.GetEntriesFast()<fDepth) return;
185
186   // pool depth reached, do mixing
187   DoMixing(pool,diele);
188
189   // increase counter for full bins
190   if (diele->fHistos) {
191     diele->fHistos->Fill("Mixing","Stats",0);
192     diele->fHistos->Fill("Mixing","CompletePools",bin);
193   }
194
195 }
196
197 //______________________________________________
198 void AliDielectronMixingHandler::DoMixing(TClonesArray &pool, AliDielectron *diele)
199 {
200   //
201   // perform the mixing
202   //
203
204   //buffer track arrays and copy them back afterwards
205   TObjArray arrTrDummy[4];
206   for (Int_t i=0; i<4; ++i) arrTrDummy[i]=diele->fTracks[i];
207
208   //buffer also global event data
209   Double_t values[AliDielectronVarManager::kNMaxValues]={0};
210   for (Int_t i=AliDielectronVarManager::kPairMax; i<AliDielectronVarManager::kNMaxValues; ++i)
211     values[i]=AliDielectronVarManager::GetValue((AliDielectronVarManager::ValueTypes)i);
212
213
214   // use event data from the first event
215   // all events are in the same mixing bin anyhow...
216   // optionally use only the event data from the first event and no tracks of it,
217   // by this you should get ride of event plane - leg correlations in the mixing
218   // but you loose the 1st event in the mixing statistics
219   AliDielectronEvent *ev0=static_cast<AliDielectronEvent*>(pool.At(0));
220   AliDielectronVarManager::SetEventData(ev0->GetEventData());
221   Int_t firstEvt=(fSkipFirstEvt ? 1 : 0);
222
223   for (Int_t i1=firstEvt; i1<pool.GetEntriesFast(); ++i1){
224     AliDielectronEvent *ev1=static_cast<AliDielectronEvent*>(pool.At(i1));
225     //use event data from the first event
226     //both events are in the same mixing bin anyhow...
227     if( !fSkipFirstEvt ) AliDielectronVarManager::SetEventData(ev1->GetEventData());
228
229     TObject *o=0x0;
230     TIter ev1P(ev1->GetTrackArrayP());
231     TIter ev1N(ev1->GetTrackArrayN());
232     
233     for (Int_t i2=i1+1; i2<pool.GetEntriesFast(); ++i2){
234       //clear arryas
235       diele->fTracks[0].Clear();
236       diele->fTracks[1].Clear();
237       diele->fTracks[2].Clear();
238       diele->fTracks[3].Clear();
239       
240       //setup track arrays
241       AliDielectronEvent *ev2=static_cast<AliDielectronEvent*>(pool.At(i2));
242       ev1P.Reset();
243       ev1N.Reset();
244       TIter ev2P(ev2->GetTrackArrayP());
245       TIter ev2N(ev2->GetTrackArrayN());
246
247       //
248       //move tracks to the same vertex (vertex of the first event), if requested
249       //
250       if (fMoveToSameVertex){
251         const Double_t *varsFirst=ev1->GetEventData();
252         const Double_t *varsMix=ev2->GetEventData();
253
254         const Double_t vFirst[3]={varsFirst[AliDielectronVarManager::kXvPrim],
255                                   varsFirst[AliDielectronVarManager::kYvPrim],
256                                   varsFirst[AliDielectronVarManager::kZvPrim]};
257
258         const Double_t vMix[3]  ={varsMix[AliDielectronVarManager::kXvPrim],
259                                   varsMix[AliDielectronVarManager::kYvPrim],
260                                   varsMix[AliDielectronVarManager::kZvPrim]};
261                                   
262         //loop over all tracks from the second event and move them to the vertex of the first
263         AliVTrack *vtrack=0x0;
264         while ( ( vtrack=(AliVTrack*)ev2P() ) ){
265           MoveToSameVertex(vtrack, vFirst, vMix);
266         }
267
268         
269         while ( ( vtrack=(AliVTrack*)ev2N() ) ){
270           MoveToSameVertex(vtrack, vFirst, vMix);
271         }
272
273         
274         ev2P.Reset();
275         ev2N.Reset();
276       }
277
278       //mixing of ev1- ev2+ (pair type4). This is common for all mixing types
279       while ( (o=ev1N()) ) diele->fTracks[1].Add(o); 
280       while ( (o=ev2P()) ) diele->fTracks[2].Add(o);
281       diele->FillPairArrays(1,2);
282       
283       if (fMixType==kAll || fMixType==kOSandLS){
284         // all 4 pair arrays will be filled
285         while ( (o=ev1P()) ) diele->fTracks[0].Add(o);
286         while ( (o=ev2N()) ) diele->fTracks[3].Add(o);
287         diele->FillPairArrays(0,2);
288         diele->FillPairArrays(1,3);
289         if (fMixType==kAll) diele->FillPairArrays(0,3);
290       }
291       
292       if (fMixType==kOSonly || fMixType==kOSandLS){
293         //use the pair type of ev1- ev1+ also for ev1+ ev1-
294         diele->fTracks[1].Clear();
295         diele->fTracks[2].Clear();
296         while ( (o=ev1P()) ) diele->fTracks[1].Add(o);
297         while ( (o=ev2N()) ) diele->fTracks[2].Add(o);
298         diele->FillPairArrays(1,2);
299       }
300
301     }
302   }
303
304   //copy back the tracks
305   for (Int_t i=0; i<4; ++i) {
306     diele->fTracks[i].Clear();
307     diele->fTracks[i]=arrTrDummy[i];
308   }
309
310   //set back global event values
311   AliDielectronVarManager::SetEventData(values);
312 }
313
314 //______________________________________________
315 Bool_t AliDielectronMixingHandler::MixRemaining(AliDielectron *diele, Int_t ipool)
316 {
317   //
318   // mix all pools even if they are incomplete
319   //
320
321   //Check if there was any processed data and it is requested to mix incomplete bins
322   if (!diele || !fMixIncomplete ) return 0;
323
324   AliDielectronVarManager::SetEvent(0x0);
325     TClonesArray *poolp=static_cast<TClonesArray*>(fArrPools.At(ipool));
326     if (!poolp || !poolp->GetEntriesFast() || !poolp->At(0)) return 0;
327     //clear the arrays before the final processing"
328     AliDebug(10,Form("Incomplete: Bin %d (%d)\n",ipool,poolp->GetEntriesFast()));
329     diele->ClearArrays();
330     DoMixing(*poolp,diele);
331
332     // increase counter for incomplete bins
333     if (diele->fHistos) {
334       //buffer event data and set event data using the first event in this pool
335       Double_t values[AliDielectronVarManager::kNMaxValues]={0};
336       for (Int_t i=AliDielectronVarManager::kPairMax; i<AliDielectronVarManager::kNMaxValues; ++i)
337         values[i]=AliDielectronVarManager::GetValue((AliDielectronVarManager::ValueTypes)i);
338       
339       AliDielectronEvent *ev1=static_cast<AliDielectronEvent*>(poolp->At(0));
340       //use event data from the first event all events are in the same mixing bin anyhow...
341       AliDielectronVarManager::SetEventData(ev1->GetEventData());
342       
343       // fill the histograms
344       diele->FillHistograms(0x0, kTRUE);
345       diele->fHistos->Fill("Mixing","Stats",1);
346       diele->fHistos->Fill("Mixing","InCompletePools",ipool);
347       diele->fHistos->Fill("Mixing","Entries_InCompletePools",poolp->GetEntriesFast());
348       
349       //set back global event values (this would mean set back to zero)
350       //AliDielectronVarManager::SetEventData(values);
351     }
352     return 1;
353 }
354
355
356 //______________________________________________
357 void AliDielectronMixingHandler::Init(const AliDielectron *diele)
358 {
359   //
360   // initialise event buffers
361   //
362
363   Int_t size=GetNumberOfBins();
364
365   AliDebug(10,Form("Creating a pool array with size %d \n",size));
366
367   if(diele->DoEventProcess()) fArrPools.Expand(size);
368
369   //add statics histogram if we have a histogram manager
370   if (diele && diele->fHistos && diele->DoEventProcess()) {
371     diele->fHistos->AddClass("Mixing");
372     diele->fHistos->UserHistogram("Mixing","Stats","Mixing Statistics;;#called bins",2,0,2);
373     TH1* h=diele->fHistos->GetHistogram("Mixing","Stats");
374     h->GetXaxis()->SetBinLabel(1,"Complete");
375     h->GetXaxis()->SetBinLabel(2,"Incomplete");
376
377     diele->fHistos->UserHistogram("Mixing","CompletePools","Mixing Statistics compete pools;bin;#fills",size,0,size);
378     diele->fHistos->UserHistogram("Mixing","InCompletePools","Mixing Statistics incomplete pools;bin;#fills",size,0,size);
379     diele->fHistos->UserHistogram("Mixing","Entries_InCompletePools","#entries in incomplete pools;entries;#fills",fDepth,0,fDepth);
380   }
381
382   TString values;
383   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i){
384     TVectorD *bins=static_cast<TVectorD*>(fAxes.At(i));
385     Int_t nRows=bins->GetNrows();
386     values+=Form("%s: ",AliDielectronVarManager::GetValueName(fEventCuts[i]));
387     for (Int_t irow=0; irow<nRows; ++irow){
388       values+=Form("%.2f, ",(*bins)[irow]);
389     }
390   }
391   
392   if (!fPID){
393     fPID=TProcessID::AddProcessID();
394   }
395
396   AliDebug(10,values.Data());
397 }
398
399 //______________________________________________
400 Int_t AliDielectronMixingHandler::GetNumberOfBins() const
401 {
402   //
403   // return the number of bins this mixing handler has
404   //
405   Int_t size=1;
406   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
407     size*=((static_cast<TVectorD*>(fAxes.At(i)))->GetNrows()-1);
408   return size;
409 }
410
411 //______________________________________________
412 Int_t AliDielectronMixingHandler::FindBin(const Double_t values[], TString *dim)
413 {
414   //
415   // bin bin in mixing stack described by 'values'
416   // if the values are outside the binning range -1 is returned
417   // if dim is non NULL debug info will be stored in the variable
418   //
419
420   if (fAxes.GetEntriesFast()==0) {
421     if (dim) (*dim)="single bin";
422     return 0;
423   }
424   if (dim) (*dim)="";
425   Int_t sizeAdd=1;
426   Int_t bin=0;
427   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i){
428     Double_t val=values[fEventCuts[i]];
429     TVectorD *bins=static_cast<TVectorD*>(fAxes.At(i));
430     Int_t nRows=bins->GetNrows();
431     if ( (val<(*bins)[0]) || (val>(*bins)[nRows-1]) ) {
432       return -1;
433     }
434
435     Int_t pos=TMath::BinarySearch(nRows,bins->GetMatrixArray(),val);
436     bin+=sizeAdd*pos;
437     if (dim) (*dim)+=Form("%s: %f (%d); ",AliDielectronVarManager::GetValueName(fEventCuts[i]),val,pos);
438     sizeAdd*=(nRows-1);
439   }
440
441   return bin;
442 }
443
444 //______________________________________________
445 void AliDielectronMixingHandler::MoveToSameVertex(AliVTrack * const vtrack, const Double_t *vFirst, const Double_t* vMix)
446 {
447   //
448   // move 'track' which belongs to the vertex information of vMix to the vertex of vFirst
449   //
450
451   static Bool_t printed=kFALSE;
452   
453   if (vtrack->IsA()==AliESDtrack::Class()){
454     AliESDtrack *track=(AliESDtrack*)vtrack;
455
456     //get track information
457     Double_t x        = track->GetX();
458     Double_t alpha    = track->GetAlpha();
459     Double_t param[5] = {0};
460     Double_t cov[15]  = {0};
461
462     for (Int_t i=0; i<5;  ++i) param[i]=track->GetParameter()[i];
463     for (Int_t i=0; i<15; ++i) cov[i]  =track->GetCovariance()[i];
464
465     //translation
466     Double_t vt[3] = {vMix[0]-vFirst[0],vMix[1]-vFirst[1],vMix[2]-vFirst[2]};
467     //rotate to the track frame
468 //     track->Global2LocalPosition(vt,track->GetAlpha());
469
470     //add to track position
471 //     x        = x       -vt[0];
472 //     param[0] = param[0]-vt[1];
473 //     param[1] = param[1]-vt[2];
474     param[1] = param[1]-vt[2];
475     
476     //set updated track information
477     track->Set(x, alpha, param, cov);
478   } else {
479     //             AliAODTrack *track=(AliAODTrack*)vtrack;
480     //             Double_t pos[3]={0};
481     //             track->GetPosition(pos);
482     //             if (pos[0]>-999.){
483       //               pos[0]=pos[0]-vMix[AliDielectronVarManager::kXvPrim]+vFirst[AliDielectronVarManager::kXvPrim];
484       //               pos[1]=pos[1]-vMix[AliDielectronVarManager::kYvPrim]+vFirst[AliDielectronVarManager::kYvPrim];
485       //               pos[2]=pos[2]-vMix[AliDielectronVarManager::kZvPrim]+vFirst[AliDielectronVarManager::kZvPrim];
486       //               track->SetPosition(pos);
487 //       AliError("Move To same vertex not yet implemented for AOD!");
488     if (!printed) {
489 //       Error("AliDielectronMixingHandler::MoveToSameVertex","Move To same vertex not yet implemented for AOD!");
490       printed=kTRUE;
491     }
492       //
493       // Not that clear how to do it. In AOD track there is fPosition, fPositionAtDCA and "TRef fProdVertex"
494       // where Xv(), Yv(), Zv() returns the position of fProdVertex
495       //
496     }
497     
498 }