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