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]={0.};
 142   AliDielectronVarManager::Fill(ev,values);
 143
 144   TString dim;
 145   Int_t bin=FindBin(values,&dim);
 146
 147   //add mixing bin to event data
 148   values[AliDielectronVarManager::kMixingBin] = bin;
 149
 150   if (bin<0){
 151     AliDebug(5,Form("Bin outside range: %s",dim.Data()));
 152     return;
 153   }
 154
 155   // get mixing pool, create it if it does not yet exist.
 156   TClonesArray *poolp=static_cast<TClonesArray*>(fArrPools.At(bin));
 157   if (!poolp){
 158     AliDebug(10,Form("New pool at %d (%s)\n",bin,dim.Data()));
 159     poolp=new(fArrPools[bin]) TClonesArray("AliDielectronEvent",fDepth);
 160   }
 161   TClonesArray &pool=*poolp;
 162
 163   AliDebug(10,Form("new event at %d: %d",bin,pool.GetEntriesFast()));
 164   AliDielectronEvent *event=new(pool[pool.GetEntriesFast()]) AliDielectronEvent();
 165   if(ev->IsA() == AliAODEvent::Class()) event->SetAOD();
 166   event->SetTracks(*diele->GetTrackArray(0), *diele->GetTrackArray(1), *diele->GetPairArray(1));
 167   event->SetEventData(values);
 168
 169   // check if pool depth is reached.
 170   if (pool.GetEntriesFast()<fDepth) return;
 171
 172   // pool depth reached, do mixing
 173   DoMixing(pool,diele);
 174
 175   // increase counter for full bins
 176   if (diele->fHistos) {
 177     diele->fHistos->Fill("Mixing","Stats",0);
 178     diele->fHistos->Fill("Mixing","CompletePools",bin);
 179   }
 180
 181   //clear the current pool
 182   pool.Clear("C");
 183 }
 184
 185 //______________________________________________
 186 void AliDielectronMixingHandler::DoMixing(TClonesArray &pool, AliDielectron *diele)
 187 {
 188   //
 189   // perform the mixing
 190   //
 191
 192   //buffer track arrays and copy them back afterwards
 193   TObjArray arrTrDummy[4];
 194   for (Int_t i=0; i<4; ++i) arrTrDummy[i]=diele->fTracks[i];
 195
 196   //buffer also global event data
 197   Double_t values[AliDielectronVarManager::kNMaxValues]={0};
 198   for (Int_t i=AliDielectronVarManager::kPairMax; i<AliDielectronVarManager::kNMaxValues; ++i)
 199     values[i]=AliDielectronVarManager::GetValue((AliDielectronVarManager::ValueTypes)i);
 200
 201   for (Int_t i1=0; i1<pool.GetEntriesFast(); ++i1){
 202     AliDielectronEvent *ev1=static_cast<AliDielectronEvent*>(pool.At(i1));
 203     //use event data from the first event
 204     //both events are in the same mixing bin anyhow...
 205     AliDielectronVarManager::SetEventData(ev1->GetEventData());
 206
 207     TObject *o=0x0;
 208     TIter ev1P(ev1->GetTrackArrayP());
 209     TIter ev1N(ev1->GetTrackArrayN());
 210
 211     for (Int_t i2=i1+1; i2<pool.GetEntriesFast(); ++i2){
 212       //clear arryas
 213       diele->fTracks[0].Clear();
 214       diele->fTracks[1].Clear();
 215       diele->fTracks[2].Clear();
 216       diele->fTracks[3].Clear();
 217
 218       //setup track arrays
 219       AliDielectronEvent *ev2=static_cast<AliDielectronEvent*>(pool.At(i2));
 220       ev1P.Reset();
 221       ev1N.Reset();
 222       TIter ev2P(ev2->GetTrackArrayP());
 223       TIter ev2N(ev2->GetTrackArrayN());
 224
 225       //
 226       //move tracks to the same vertex (vertex of the first event), if requested
 227       //
 228       if (fMoveToSameVertex){
 229         const Double_t *varsFirst=ev1->GetEventData();
 230         const Double_t *varsMix=ev2->GetEventData();
 231
 232         const Double_t vFirst[3]={varsFirst[AliDielectronVarManager::kXvPrim],
 233                                   varsFirst[AliDielectronVarManager::kYvPrim],
 234                                   varsFirst[AliDielectronVarManager::kZvPrim]};
 235
 236         const Double_t vMix[3]  ={varsMix[AliDielectronVarManager::kXvPrim],
 237                                   varsMix[AliDielectronVarManager::kYvPrim],
 238                                   varsMix[AliDielectronVarManager::kZvPrim]};
 239
 240         //loop over all tracks from the second event and move them to the vertex of the first
 241         AliVTrack *vtrack=0x0;
 242         while ( ( vtrack=(AliVTrack*)ev2P() ) ){
 243           MoveToSameVertex(vtrack, vFirst, vMix);
 244         }
 245
 246
 247         while ( ( vtrack=(AliVTrack*)ev2N() ) ){
 248           MoveToSameVertex(vtrack, vFirst, vMix);
 249         }
 250
 251
 252         ev2P.Reset();
 253         ev2N.Reset();
 254       }
 255
 256       // set the vertex to the track in case (AOD only)
 257       if ( ev1->IsAOD() && ev2->IsAOD() ) {
 258         Int_t idx = 0;
 259
 260         while ( (o=ev1N()) ) {
 261           idx = ev1->GetTrackArrayN()->IndexOf(o);
 262           ((AliAODTrack *)o)->SetProdVertex(ev1->GetVertexArrayN()->At(idx));
 263         }
 264         while ( (o=ev1P()) ) {
 265           idx = ev1->GetTrackArrayP()->IndexOf(o);
 266           ((AliAODTrack *)o)->SetProdVertex(ev1->GetVertexArrayP()->At(idx));
 267         }
 268         while ( (o=ev2N()) ) {
 269           idx = ev2->GetTrackArrayN()->IndexOf(o);
 270           ((AliAODTrack *)o)->SetProdVertex(ev2->GetVertexArrayN()->At(idx));
 271         }
 272         while ( (o=ev2P()) ) {
 273           idx = ev2->GetTrackArrayP()->IndexOf(o);
 274           ((AliAODTrack *)o)->SetProdVertex(ev2->GetVertexArrayP()->At(idx));
 275         }
 276
 277         // reset the iterators
 278         ev1N.Reset();
 279         ev1P.Reset();
 280         ev2N.Reset();
 281         ev2P.Reset();
 282       }
 283
 284       //mixing of ev1- ev2+ (pair type4). This is common for all mixing types
 285       while ( (o=ev1N()) ) diele->fTracks[1].Add(o);
 286       while ( (o=ev2P()) ) diele->fTracks[2].Add(o);
 287       diele->FillPairArrays(1,2);
 288
 289       if (fMixType==kAll || fMixType==kOSandLS){
 290         // all 4 pair arrays will be filled
 291         while ( (o=ev1P()) ) diele->fTracks[0].Add(o);
 292         while ( (o=ev2N()) ) diele->fTracks[3].Add(o);
 293         diele->FillPairArrays(0,2);
 294         diele->FillPairArrays(1,3);
 295         if (fMixType==kAll) diele->FillPairArrays(0,3);
 296       }
 297
 298       if (fMixType==kOSonly || fMixType==kOSandLS){
 299         //use the pair type of ev1- ev1+ also for ev1+ ev1-
 300         diele->fTracks[1].Clear();
 301         diele->fTracks[2].Clear();
 302         while ( (o=ev1P()) ) diele->fTracks[1].Add(o);
 303         while ( (o=ev2N()) ) diele->fTracks[2].Add(o);
 304         diele->FillPairArrays(1,2);
 305       }
 306
 307     }
 308   }
 309
 310   //copy back the tracks
 311   for (Int_t i=0; i<4; ++i) {
 312     diele->fTracks[i].Clear();
 313     diele->fTracks[i]=arrTrDummy[i];
 314   }
 315
 316   //set back global event values
 317   AliDielectronVarManager::SetEventData(values);
 318 }
 319
 320 //______________________________________________
 321 void AliDielectronMixingHandler::MixRemaining(AliDielectron *diele)
 322 {
 323   //
 324   // mix all pools even if they are incomplete
 325   //
 326
 327   //Check if there was any processed data and it is requested to mix incomplete bins
 328   if (!diele || !diele->PairArray(0) || !fMixIncomplete ) return;
 329
 330   AliDielectronVarManager::SetEvent(0x0);
 331   for (Int_t ipool=0; ipool<fArrPools.GetSize(); ++ipool){
 332     TClonesArray *poolp=static_cast<TClonesArray*>(fArrPools.At(ipool));
 333     if (!poolp || !poolp->GetEntriesFast() || !poolp->At(0)) continue;
 334     //clear the arrays before the final processing"
 335     AliDebug(10,Form("Incomplete: Bin %d (%d)\n",ipool,poolp->GetEntriesFast()));
 336     diele->ClearArrays();
 337     DoMixing(*poolp,diele);
 338
 339     // increase counter for incomplete bins
 340     if (diele->fHistos) {
 341       //buffer event data and set event data using the first event in this pool
 342       Double_t values[AliDielectronVarManager::kNMaxValues]={0};
 343       for (Int_t i=AliDielectronVarManager::kPairMax; i<AliDielectronVarManager::kNMaxValues; ++i)
 344         values[i]=AliDielectronVarManager::GetValue((AliDielectronVarManager::ValueTypes)i);
 345
 346       AliDielectronEvent *ev1=static_cast<AliDielectronEvent*>(poolp->At(0));
 347       //use event data from the first event all events are in the same mixing bin anyhow...
 348       AliDielectronVarManager::SetEventData(ev1->GetEventData());
 349
 350       // fill the histograms
 351       diele->FillHistograms(0x0, kTRUE);
 352       diele->fHistos->Fill("Mixing","Stats",1);
 353       diele->fHistos->Fill("Mixing","InCompletePools",ipool);
 354       diele->fHistos->Fill("Mixing","Entries_InCompletePools",poolp->GetEntriesFast());
 355
 356       //set back global event values
 357       AliDielectronVarManager::SetEventData(values);
 358     }
 359
 360   }
 361 }
 362
 363
 364 //______________________________________________
 365 void AliDielectronMixingHandler::Init(const AliDielectron *diele)
 366 {
 367   //
 368   // initialise event buffers
 369   //
 370
 371   Int_t size=GetNumberOfBins();
 372
 373   AliDebug(10,Form("Creating a pool array with size %d \n",size));
 374
 375   fArrPools.Expand(size);
 376
 377   //add statics histogram if we have a histogram manager
 378   if (diele && diele->fHistos) {
 379     diele->fHistos->AddClass("Mixing");
 380     diele->fHistos->UserHistogram("Mixing","Stats","Mixing Statistics;;#called bins",2,0,2);
 381     TH1* h=diele->fHistos->GetHistogram("Mixing","Stats");
 382     h->GetXaxis()->SetBinLabel(1,"Complete");
 383     h->GetXaxis()->SetBinLabel(2,"Incomplete");
 384
 385     diele->fHistos->UserHistogram("Mixing","CompletePools","Mixing Statistics compete pools;bin;#fills",size,0,size);
 386     diele->fHistos->UserHistogram("Mixing","InCompletePools","Mixing Statistics incomplete pools;bin;#fills",size,0,size);
 387     diele->fHistos->UserHistogram("Mixing","Entries_InCompletePools","#entries in incomplete pools;entries;#fills",fDepth,0,fDepth);
 388   }
 389
 390   TString values;
 391   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i){
 392     TVectorD *bins=static_cast<TVectorD*>(fAxes.At(i));
 393     Int_t nRows=bins->GetNrows();
 394     values+=Form("%s: ",AliDielectronVarManager::GetValueName(fEventCuts[i]));
 395     for (Int_t irow=0; irow<nRows; ++irow){
 396       values+=Form("%.2f, ",(*bins)[irow]);
 397     }
 398   }
 399
 400   AliDebug(10,values.Data());
 401 }
 402
 403 //______________________________________________
 404 Int_t AliDielectronMixingHandler::GetNumberOfBins() const
 405 {
 406   //
 407   // return the number of bins this mixing handler has
 408   //
 409   Int_t size=1;
 410   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
 411     size*=((static_cast<TVectorD*>(fAxes.At(i)))->GetNrows()-1);
 412   return size;
 413 }
 414
 415 //______________________________________________
 416 Int_t AliDielectronMixingHandler::FindBin(const Double_t values[], TString *dim)
 417 {
 418   //
 419   // bin bin in mixing stack described by 'values'
 420   // if the values are outside the binning range -1 is returned
 421   // if dim is non NULL debug info will be stored in the variable
 422   //
 423
 424   if (fAxes.GetEntriesFast()==0) {
 425     if (dim) (*dim)="single bin";
 426     return 0;
 427   }
 428   if (dim) (*dim)="";
 429   Int_t sizeAdd=1;
 430   Int_t bin=0;
 431   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i){
 432     Double_t val=values[fEventCuts[i]];
 433     TVectorD *bins=static_cast<TVectorD*>(fAxes.At(i));
 434     Int_t nRows=bins->GetNrows();
 435     if ( (val<(*bins)[0]) || (val>(*bins)[nRows-1]) ) {
 436       return -1;
 437     }
 438
 439     Int_t pos=TMath::BinarySearch(nRows,bins->GetMatrixArray(),val);
 440     bin+=sizeAdd*pos;
 441     if (dim) (*dim)+=Form("%s: %f (%d); ",AliDielectronVarManager::GetValueName(fEventCuts[i]),val,pos);
 442     sizeAdd*=(nRows-1);
 443   }
 444
 445   return bin;
 446 }
 447
 448 //______________________________________________
 449 void AliDielectronMixingHandler::MoveToSameVertex(AliVTrack * const vtrack, const Double_t *vFirst, const Double_t* vMix)
 450 {
 451   //
 452   // move 'track' which belongs to the vertex information of vMix to the vertex of vFirst
 453   //
 454
 455   static Bool_t printed=kFALSE;
 456
 457   if (vtrack->IsA()==AliESDtrack::Class()){
 458     AliESDtrack *track=(AliESDtrack*)vtrack;
 459
 460     //get track information
 461     Double_t x        = track->GetX();
 462     Double_t alpha    = track->GetAlpha();
 463     Double_t param[5] = {0};
 464     Double_t cov[15]  = {0};
 465
 466     for (Int_t i=0; i<5;  ++i) param[i]=track->GetParameter()[i];
 467     for (Int_t i=0; i<15; ++i) cov[i]  =track->GetCovariance()[i];
 468
 469     //translation
 470     Double_t vt[3] = {vMix[0]-vFirst[0],vMix[1]-vFirst[1],vMix[2]-vFirst[2]};
 471     //rotate to the track frame
 472 //     track->Global2LocalPosition(vt,track->GetAlpha());
 473
 474     //add to track position
 475 //     x        = x       -vt[0];
 476 //     param[0] = param[0]-vt[1];
 477 //     param[1] = param[1]-vt[2];
 478     param[1] = param[1]-vt[2];
 479
 480     //set updated track information
 481     track->Set(x, alpha, param, cov);
 482   } else {
 483     //             AliAODTrack *track=(AliAODTrack*)vtrack;
 484     //             Double_t pos[3]={0};
 485     //             track->GetPosition(pos);
 486     //             if (pos[0]>-999.){
 487       //               pos[0]=pos[0]-vMix[AliDielectronVarManager::kXvPrim]+vFirst[AliDielectronVarManager::kXvPrim];
 488       //               pos[1]=pos[1]-vMix[AliDielectronVarManager::kYvPrim]+vFirst[AliDielectronVarManager::kYvPrim];
 489       //               pos[2]=pos[2]-vMix[AliDielectronVarManager::kZvPrim]+vFirst[AliDielectronVarManager::kZvPrim];
 490       //               track->SetPosition(pos);
 491 //       AliError("Move To same vertex not yet implemented for AOD!");
 492     if (!printed) {
 493 //       Error("AliDielectronMixingHandler::MoveToSameVertex","Move To same vertex not yet implemented for AOD!");
 494       printed=kTRUE;
 495     }
 496       //
 497       // Not that clear how to do it. In AOD track there is fPosition, fPositionAtDCA and "TRef fProdVertex"
 498       // where Xv(), Yv(), Zv() returns the position of fProdVertex
 499       //
 500     }
 501
 502 }