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