PWG1/background/AliAnalysisTaskBGvsTime.cxx

   1 // AliAnalysisTaskBGvsTime
   2
   3 // This task computes various histos as a function of time, using a
   4 // time stamp made from orbit, bx number and period. It is used for
   5 // studies of BG vs time and of luminosity.
   6
   7 // The histograms are booked dinamically when they are called the
   8 // first time. This simplifies the handling of different trigger
   9 // classes (you don't need to know them in advance). The merging is
  10 // then complicated, but it is handled via the external class
  11 // AliHistoListWrapper.
  12
  13 // Different sets of cuts can be used, and the binning of vs time
  14 // histos can be defined dynamically (see setters)
  15
  16 // Can also process MC, in which case the efficiency x acceptance for
  17 // tracks used in the rate (luminosity) histos is also computed, and
  18 // put in a histo in which each bin corresponds to a given rate histo
  19 // (eta < 0.8 and pt > 0.5 or pt > 1.0 at the time I'm writing this).
  20
  21 // The rates are corrected for the dead time, computed using trigger
  22 // scalers. WARNING: the dead time computation will NOT WORK on CAF
  23 // (it assumes consecutive events).
  24
  25 // Author: Michele Floris, CERN
  26
  27 #include "AliAnalysisTaskBGvsTime.h"
  28 #include "AliESDInputHandler.h"
  29 #include "TString.h"
  30 #include "AliVParticle.h"
  31 #include "AliESDInputHandlerRP.h"
  32 #include "TTree.h"
  33 #include "TClonesArray.h"
  34 #include "TBranch.h"
  35 //#include "AliITSRecPoint.h"
  36 #include "AliMultiplicity.h"
  37 #include "AliLog.h"
  38 #include "AliAnalysisManager.h"
  39 #include "AliBackgroundSelection.h"
  40 #include <iostream>
  41 #include "TFile.h"
  42 #include "TCanvas.h"
  43 #include "AliHistoListWrapper.h"
  44 #include "AliTriggerAnalysis.h"
  45 #include "TMath.h"
  46 #include "AliPhysicsSelection.h"
  47 #include "AliBackgroundSelection.h"
  48 #include "AliESDtrackCuts.h"
  49 #include "AliStack.h"
  50 #include "AliMCEvent.h"
  51 #include "TDatime.h"
  52
  53 using namespace std;
  54
  55 ClassImp(AliAnalysisTaskBGvsTime)
  56
  57
  58
  59 AliAnalysisTaskBGvsTime::AliAnalysisTaskBGvsTime()
  60   : AliAnalysisTaskSE("TaskBGvsTime"),
  61     fESD(0),fListHisto(0),fListWrapper(0),fStartTime(0),fEndTime(0),
  62     fNMultBins(0), fMultBins(0),fFirstTimeStamp(0), fLastTimeStamp(0), fBinWidth(10), fNoCuts(0), fPhysicsSelection(0), fUsePhysicsSelection(0),
  63     fUseZeroBin(0), fIsMC(0), fSkipV0(0), fSkipZeroBin(0), fUseBunchInt(0), fHistoTimeStampVsUTC(0), fHistoTimeStampDiffUTC(0)
  64
  65 {
  66   // constructor
  67
  68   DefineOutput(1, AliHistoListWrapper::Class());
  69   DefineOutput(2, AliPhysicsSelection::Class());
  70   //  DefineOutput(2, TH1I::Class());
  71
  72 }
  73 AliAnalysisTaskBGvsTime::AliAnalysisTaskBGvsTime(const char * name)
  74   : AliAnalysisTaskSE(name),fESD (0),fListHisto(0),fListWrapper(0),fStartTime(0),fEndTime(0),
  75     fNMultBins(0),fMultBins(0),fFirstTimeStamp(0), fLastTimeStamp(0), fBinWidth(10), fNoCuts(0), fPhysicsSelection(0), fUsePhysicsSelection(0),
  76     fUseZeroBin(0), fIsMC(0), fSkipV0(0), fSkipZeroBin(0), fUseBunchInt(0), fHistoTimeStampVsUTC(0), fHistoTimeStampDiffUTC(0)
  77 {
  78   //
  79   // Standard constructur which should be used
  80   //
  81
  82   DefineOutput(1, AliHistoListWrapper::Class());
  83   DefineOutput(2, AliPhysicsSelection::Class());
  84   //  DefineOutput(2, TH1I::Class());
  85
  86 }
  87
  88 AliAnalysisTaskBGvsTime::AliAnalysisTaskBGvsTime(const AliAnalysisTaskBGvsTime& obj) :
  89   AliAnalysisTaskSE(obj) ,fESD (0),fListHisto(0),fListWrapper(0),fStartTime(0),fEndTime(0),
  90   fNMultBins(0),fMultBins(0),fFirstTimeStamp(0), fLastTimeStamp(0), fBinWidth(10), fNoCuts(0), fPhysicsSelection(0), fUsePhysicsSelection(0),
  91   fUseZeroBin(0), fIsMC(0), fSkipV0(0), fSkipZeroBin(0), fUseBunchInt(0), fHistoTimeStampVsUTC(0), fHistoTimeStampDiffUTC(0)
  92 {
  93   //copy ctor
  94   fESD = obj.fESD ;
  95   fListHisto= obj.fListHisto;
  96   fListWrapper= obj.fListWrapper;
  97   fStartTime= obj.fStartTime;
  98   fEndTime= obj.fEndTime;
  99
 100   fNMultBins= obj.fNMultBins;
 101   fMultBins= obj.fMultBins;
 102   fFirstTimeStamp= obj.fFirstTimeStamp;
 103   fLastTimeStamp= obj.fLastTimeStamp;
 104   fBinWidth= obj.fBinWidth;
 105   fNoCuts= obj.fNoCuts;
 106   fPhysicsSelection= obj.fPhysicsSelection;
 107   fUsePhysicsSelection= obj.fUsePhysicsSelection;
 108   fUseZeroBin= obj.fUseZeroBin;
 109   fIsMC= obj.fIsMC;
 110   fSkipV0= obj.fSkipV0;
 111   fSkipZeroBin= obj.fSkipZeroBin;
 112   fUseBunchInt= obj.fUseBunchInt;
 113   fHistoTimeStampVsUTC = obj.fHistoTimeStampVsUTC;
 114   fHistoTimeStampDiffUTC = obj.fHistoTimeStampDiffUTC;
 115
 116 }
 117
 118 AliAnalysisTaskBGvsTime::~AliAnalysisTaskBGvsTime(){
 119
 120   // destructor
 121   if(fListWrapper) {
 122     delete fListWrapper;
 123     fListWrapper = 0;
 124   }
 125   if(fMultBins) {
 126     delete fMultBins;
 127     fMultBins = 0;
 128   }
 129   if(fPhysicsSelection) {
 130     delete fPhysicsSelection;
 131     fPhysicsSelection = 0;
 132   }
 133   // Histos should not be destroyed: fListWrapper is owner!
 134 }
 135 void AliAnalysisTaskBGvsTime::UserCreateOutputObjects()
 136 {
 137   // Called once
 138   fListWrapper = new AliHistoListWrapper("histoListWrapper","histoListWrapper");
 139   fListHisto = fListWrapper->GetList();
 140   Float_t lenght =  fEndTime - fStartTime;
 141   Int_t nBins = TMath::FloorNint(lenght/fBinWidth)+1;
 142   Int_t end   = nBins*fBinWidth;
 143   // cover 20 days in UTC. I'll set offset by hand below.
 144   Int_t lenghtUTC = 20*3600*24;
 145   Int_t nBinsUTC = TMath::FloorNint(Double_t(lenghtUTC)/fBinWidth)+1;
 146   Int_t endUTC   = nBinsUTC*fBinWidth;
 147   if (!fIsMC) {
 148     Bool_t oldStatus = TH1::AddDirectoryStatus();
 149     TH1::AddDirectory(kFALSE);
 150     fHistoTimeStampVsUTC = new TH2F ("fHistoTimeStampVsUTC", "fHistoTimeStampVsUTC", nBins, -0.5, end-0.5,  nBinsUTC, -0.5, endUTC-0.5);
 151     fListHisto->Add(fHistoTimeStampVsUTC);
 152     fHistoTimeStampDiffUTC = new TH1F("fHistoTimeStampDiffUTC", "fHistoTimeStampDiffUTC", nBinsUTC, -0.5, endUTC-0.5);
 153     fListHisto->Add(fHistoTimeStampDiffUTC);
 154     TH1::AddDirectory(oldStatus);
 155   }
 156
 157   fPhysicsSelection = new AliPhysicsSelection();
 158   fPhysicsSelection->SetUseBXNumbers();
 159   if(fIsMC) fPhysicsSelection->SetAnalyzeMC();
 160   fPhysicsSelection->SetBin0Callback(this->GetName());
 161
 162   fPhysicsSelection->SetSkipTriggerClassSelection();//
 163
 164   AliBackgroundSelection * bg = new AliBackgroundSelection();
 165   bg->SetDeltaPhiCut(10);
 166   fPhysicsSelection->AddBackgroundIdentification(bg);
 167
 168   if (fSkipV0) {
 169     fPhysicsSelection->SetSkipV0();
 170   }
 171
 172 }
 173
 174
 175 void AliAnalysisTaskBGvsTime::UserExec(Option_t *)
 176 {
 177
 178   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 179   PostData(1,fListWrapper);
 180   PostData(2,fPhysicsSelection);
 181
 182
 183   const float etaCut = 0.8;
 184
 185
 186   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 187   if(!fESD) {
 188     AliFatal("Cannot get ESD");
 189   }
 190   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
 191     AliFatal("Not processing ESDs");
 192   }
 193
 194   if (fUseBunchInt) {
 195     fPhysicsSelection->SetComputeBG();
 196   }
 197
 198   // Get V0 flags and trigger flags
 199   static AliTriggerAnalysis * triggerAnalysis = new AliTriggerAnalysis();
 200   Bool_t v0A   = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0A);
 201   Bool_t v0C   = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0C);
 202   Bool_t v0ABG = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0ABG);
 203   Bool_t v0CBG = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0CBG);
 204   Bool_t mb1Offline = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kMB1);
 205   Bool_t isInV0 = !v0ABG && !v0CBG && ((v0A && !v0C) || (v0C && !v0A)); // try to select beam gas in CINT1A/C events: require one v0 interaction (but not the other) and not BG hit in v0. This should select events produced in between the 2 v0s and boosted forward.
 206
 207
 208   // If it is MC: fill histo of generated events for efficiency calculations
 209   if (fIsMC) {
 210     Bool_t atLeastPt1MC  = kFALSE;
 211     Bool_t atLeastPt05MC = kFALSE;
 212
 213     if (!fMCEvent) {
 214       AliError("No MC info found");
 215     } else {
 216
 217       //loop on the MC event
 218       for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 219         AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 220
 221         // We don't care about neutrals and non-physical primaries
 222         if(mcPart->Charge() == 0) continue;
 223         if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;
 224
 225         // Kinematic cuts:
 226         if (TMath::Abs(mcPart->Eta()) < etaCut) {
 227           if (mcPart->Pt() > 0.5) atLeastPt05MC  = kTRUE;
 228           if (mcPart->Pt() > 1.0) atLeastPt1MC   = kTRUE;
 229         }
 230         if (atLeastPt1MC && atLeastPt05MC) break; // no need to look for other tracks
 231       }
 232       if (atLeastPt1MC ) {
 233         GetEfficiencyHisto(kEffStepGen)->Fill(kEffPt1);
 234         if(mb1Offline) GetEfficiencyHisto(kEffStepTrig)->Fill(kEffPt1);
 235       }
 236       if (atLeastPt05MC) {
 237         GetEfficiencyHisto(kEffStepGen)->Fill(kEffPt05 );
 238         if(mb1Offline) GetEfficiencyHisto(kEffStepTrig)->Fill(kEffPt05);
 239       }
 240     }
 241   }
 242
 243
 244
 245
 246   // CUTS
 247   const AliMultiplicity* mult = fESD->GetMultiplicity();
 248   if (!mult){
 249     AliFatal("No multiplicity object"); // TODO: Should this be fatal?
 250   }
 251   // get number of SPD clusters
 252   Int_t spdClusters = 0;
 253   for(Int_t ilayer = 0; ilayer < 2; ilayer++){
 254     spdClusters += mult->GetNumberOfITSClusters(ilayer);
 255   }
 256
 257   //  Bool_t isZeroBin = kFALSE;
 258   Bool_t isZeroBin = IsEventInBinZero();
 259   if(fUseZeroBin  && !isZeroBin) return;
 260   if(fSkipZeroBin && isZeroBin ) return;
 261
 262   Bool_t physelDecision = fPhysicsSelection->IsCollisionCandidate(fESD);
 263
 264   if (fUsePhysicsSelection) {
 265     if(!physelDecision) return;
 266   }
 267   else if (!fNoCuts) {if (spdClusters < 2) return;}// At least 2 clusters
 268   //  else {AliInfo("No Cuts");}
 269
 270
 271
 272   // get time stamp
 273   Long64_t timeStampBX = 0;
 274   timeStampBX = fESD->GetBunchCrossNumber();
 275   timeStampBX += (Long64_t) 3564 * (fESD->GetOrbitNumber() + fESD->GetPeriodNumber() * 16777216);
 276   timeStampBX = (Long64_t) (25e-9 * timeStampBX);
 277   if (fFirstTimeStamp == 0) {
 278     fFirstTimeStamp = timeStampBX;
 279     fLastTimeStamp  = timeStampBX;
 280   }
 281   if (timeStampBX < fFirstTimeStamp) {
 282     AliError("Time stamp not monothonic!");
 283     fFirstTimeStamp = timeStampBX;
 284   }
 285   if (timeStampBX > fLastTimeStamp) {
 286     fLastTimeStamp = timeStampBX;
 287   }
 288   timeStampBX -= fStartTime;
 289
 290   Long64_t timeStamp = timeStampBX;
 291
 292   Long64_t timeStampGDC = fESD->GetTimeStamp();
 293
 294   static TDatime timeOffsetDT(2009,12,5,0,0,0);
 295   static Long64_t timeOffset = timeOffsetDT.Convert();
 296
 297
 298
 299   // Get trigger scalers for dead time calculation (only data)
 300   // Only CINT1B (at least for the time being)
 301   AliESDHeader* esdheader = (AliESDHeader*)fESD->GetHeader();
 302   static ULong64_t L0 = 0;
 303   static ULong64_t L2 = 0;
 304   if (fESD->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")&&!fIsMC) {
 305     AliTriggerScalersRecordESD* scalrecord = (AliTriggerScalersRecordESD*)esdheader->GetTriggerScalersRecord();
 306     const AliTriggerScalersESD* scalers = scalrecord->GetTriggerScalersForClass(2); //2 is the cint1b class index in the trigger mask
 307     L0 = scalers->GetLOCB(); //L0 before any vetos
 308     L2 = scalers->GetL2CA(); //L2 after vetos
 309   }
 310
 311
 312   // loop over trigger classes in the event
 313   TObjArray * tokens = 0;
 314   if(fIsMC) {
 315     // in case of montecarlo I override the trigger class to CINT1B for latter compatibility
 316     tokens = new TObjArray;
 317     tokens->SetOwner();
 318     tokens->Add(new TObjString("CINT1B-ABCE-NOPF-ALL"));
 319   }
 320   else {
 321     TString trgClasses = fESD->GetFiredTriggerClasses();
 322     tokens = trgClasses.Tokenize(" ");
 323   }
 324   TIter iter(tokens);
 325
 326   while(TObjString * tok = (TObjString*) iter.Next()){
 327     // clean up trigger name
 328     TString trg = tok->GetString();
 329     trg.Strip(TString::kTrailing, ' ');
 330     trg.Strip(TString::kLeading, ' ');
 331     // print selected events in !CIN1B trigs:
 332     //    if(!trg.Contains("CINT1B"))
 333 //       Printf("File: %s, IEV: %d, TRG: %s, Orbit: 0x%x, Period: %d, BC: %d\n",
 334 //                                     ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(),
 335 //                                     trg.Data(),
 336 //                                     fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
 337
 338
 339
 340     // Fill histos
 341     GetVsTimeHistoAll(trg.Data())->Fill(timeStamp);
 342     GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"ALL")->Fill(timeStamp);
 343     GetVsTimeHisto(trg.Data(),-1,         fESD->GetBunchCrossNumber(),"ALL")->Fill(timeStamp);
 344     GetVsTimeHisto(trg.Data(),spdClusters,-1,                         "ALL")->Fill(timeStamp);
 345     GetVsTimeHisto(trg.Data(),-1         ,-1,                         "ALL")->Fill(timeStamp);
 346     if (isInV0)     {
 347       GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"inV0")->Fill(timeStamp);
 348       GetVsTimeHisto(trg.Data(),-1         ,fESD->GetBunchCrossNumber(),"inV0")->Fill(timeStamp);
 349       GetVsTimeHisto(trg.Data(),-1         ,-1,                         "inV0")->Fill(timeStamp);
 350     }
 351
 352     // In order to compute mean multiplicity of spd clusters in a time
 353     // bin, we integrate the multiplicity in that bin, and then we
 354     // divide for the number of events in that bin in terminate.
 355
 356     // Is the error computed correctly if we fill with a weight ? Looping to be sure...
 357     for (Int_t iclus = 0; iclus < spdClusters; iclus ++)  {
 358       GetVsTimeHisto((TString("hMultSPDvsTime_")+trg).Data())->Fill(timeStamp);
 359       GetVsTimeHisto(trg.Data(),-1,         fESD->GetBunchCrossNumber(),"ALL","MultSPDvsTime")->Fill(timeStamp);
 360       GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"ALL","MultSPDvsTime")->Fill(timeStamp);
 361       GetVsTimeHisto(trg.Data(),spdClusters,-1,                         "ALL","MultSPDvsTime")->Fill(timeStamp);
 362       GetVsTimeHisto(trg.Data(),-1,-1,                                  "ALL","MultSPDvsTime")->Fill(timeStamp);
 363       if (isInV0)     {
 364         GetVsTimeHisto(trg.Data(),-1,         fESD->GetBunchCrossNumber(),"inV0","MultSPDvsTime")->Fill(timeStamp);
 365         GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"inV0","MultSPDvsTime")->Fill(timeStamp);
 366         GetVsTimeHisto(trg.Data(),-1,-1,                                  "inV0","MultSPDvsTime")->Fill(timeStamp);
 367       }
 368
 369     }
 370
 371     // Keep cluster distribution for reference
 372     GetDistributionHisto(trg.Data(),kDistSPDMult)->Fill(spdClusters);
 373     if(isInV0) GetDistributionHisto(trg.Data(),kDistSPDMult,"_inV0")->Fill(spdClusters);
 374
 375     // Distribution of hits per its layer:
 376     for(Int_t ilayer = 0; ilayer < 6; ilayer++){
 377       GetDistributionHisto(trg.Data(),kDistClsITSLayer)->Fill(ilayer,mult->GetNumberOfITSClusters(ilayer));// fill weighting with the number of CLS
 378     }
 379
 380
 381     // Dead time vs time stamp & time offset
 382     // WARNING THIS WON'T WORK ON CAF, NOR GRID: REDO IT WITH A NEW MERGEABLE OBJECT
 383     if (trg=="CINT1B-ABCE-NOPF-ALL"&&!fIsMC) {
 384
 385       // Fill time difference histos (only for CINt1B)
 386       if (!fIsMC) {
 387         fHistoTimeStampVsUTC->Fill(timeStampBX,timeStampGDC-timeOffset);
 388         fHistoTimeStampDiffUTC->Fill(timeStampGDC-timeOffset-timeStampBX);
 389       }
 390
 391       static ULong64_t oldL0   = 0;  // L0 counts at the beginning of this bin
 392       static ULong64_t oldL2   = 0;  // L2 counts at the beginning of this bin
 393       static ULong64_t prevL0   = 0; // L0 counts in the previous event
 394       static ULong64_t prevL2   = 0; // L2 counts in the previous event
 395
 396       static ULong64_t oldTime = 0;  // time stamp at the beginning of this bin
 397       //      static ULong64_t previousTime; // timestamp in the previous event
 398
 399       static Int_t prevbin = -1; // bin in the previous event
 400
 401       Int_t bin = GetDeadTimeHisto(trg.Data())->FindBin(timeStamp);
 402
 403       if (prevbin == -1) { // first event
 404         prevbin = bin;
 405         oldL0  = L0;
 406         oldL2  = L2;
 407         oldTime = timeStamp;
 408
 409       } else if (prevbin != bin) {
 410         // New bin: let's fill the previous one
 411         Double_t dL0 = Double_t(prevL0 - oldL0);
 412         Double_t dL2 = Double_t(prevL2 - oldL2);
 413
 414         //      Double_t deadtime  =  Double_t(1 - dL2/dL0); // interested in relative fraction of dead time
 415         Double_t deadtime  =  Double_t(dL2/dL0); // interested in relative fraction of dead time
 416         Double_t edeadtime = TMath::Sqrt(dL2*(dL0-dL2)/dL0/dL0/dL0); // Binomial error
 417
 418 //      cout << "DEADTIME " << endl;
 419 //      cout << L0 << " " << dL0 << " " << oldL0 << " " << prevL0 << endl;
 420 //      cout << L2 << " " << dL2 << " " << oldL2 << " " << prevL2 << endl;
 421 //      cout << deadtime << endl;
 422
 423
 424         GetDeadTimeHisto(trg.Data())->SetBinContent(prevbin, deadtime );
 425         GetDeadTimeHisto(trg.Data())->SetBinError  (prevbin, edeadtime);
 426
 427         oldL0  = L0;
 428         oldL2  = L2;
 429         oldTime = timeStamp;
 430
 431       }
 432       prevbin = bin;
 433       prevL0  = L0;
 434       prevL2  = L2;
 435
 436     }
 437
 438     // TPC track multiplicity vs time
 439
 440     // Selection by andrea.dainese@pd.infn.it
 441     // 15.03.2010
 442
 443
 444     // Primary vertex
 445     Bool_t badVertex = kFALSE;
 446     const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 447     if(vertex->GetNContributors()<1) {
 448       // SPD vertex
 449       vertex = fESD->GetPrimaryVertexSPD();
 450       if(vertex->GetNContributors()<1) {
 451         badVertex = kTRUE;
 452       }
 453     }
 454
 455     // Fill vertex distribution
 456     if ( ((vertex->IsFromVertexerZ() && vertex->GetDispersion()<=0.02) || !vertex->IsFromVertexerZ()) && !badVertex ) {
 457       GetDistributionHisto(trg.Data(),kDistVertex)->Fill(vertex->GetZ());
 458       if (vertex->IsFromVertexerZ()) GetDistributionHisto(trg.Data(),kDistVertexZ) ->Fill(vertex->GetZ());
 459       else                           GetDistributionHisto(trg.Data(),kDistVertex3D)->Fill(vertex->GetZ());
 460
 461       if(isInV0) GetDistributionHisto(trg.Data(),kDistVertex,"_inV0")->Fill(vertex->GetZ());
 462     }
 463
 464     // apply a cut |zVertex| < CUT, if needed
 465
 466     // Track cuts (except d0 cut)
 467     //------- TPC track selection --------
 468     // Selection by andrea.dainese@pd.infn.it
 469     // 15.03.2010
 470
 471     Bool_t selectPrimaries=kTRUE;
 472     static AliESDtrackCuts* esdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
 473
 474     // loop on tracks
 475     Int_t ntracks = fESD->GetNumberOfTracks();
 476     if(badVertex) ntracks = -1; // skip loop if the vertex is bad
 477
 478     // flags for luminosity histos
 479     Bool_t atLeastPt1  = kFALSE;
 480     Bool_t atLeastPt05 = kFALSE;
 481
 482     for (Int_t iTrack = 0; iTrack<ntracks; iTrack++) {
 483       AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTrack));
 484       // for each track
 485
 486       // track quality cuts
 487       if(!esdtrackCutsITSTPC->AcceptTrack(track)) continue;
 488
 489       // bring it to the primary vertex and compute impact parameters
 490       if(!track->RelateToVertex(vertex,fESD->GetMagneticField(),kVeryBig)) continue; // this is already done in AliReconstruction...
 491
 492       // track-to-vertex cut (see below)
 493       if(!SelectOnImpPar(track)) continue;
 494
 495       // Fill histos (TPC Multiplicity vs TIME)
 496       GetVsTimeHisto(trg.Data(),-1,         fESD->GetBunchCrossNumber(),"ALL","MultTPCvsTime")->Fill(timeStamp);
 497       GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"ALL","MultTPCvsTime")->Fill(timeStamp);
 498       GetVsTimeHisto(trg.Data(),spdClusters,-1,                         "ALL","MultTPCvsTime")->Fill(timeStamp);
 499       GetVsTimeHisto(trg.Data(),-1,-1,                                  "ALL","MultTPCvsTime")->Fill(timeStamp);
 500       if ((v0A || v0C) && !(v0ABG || v0CBG))     {
 501         GetVsTimeHisto(trg.Data(),-1,         fESD->GetBunchCrossNumber(),"inV0","MultTPCvsTime")->Fill(timeStamp);
 502         GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"inV0","MultTPCvsTime")->Fill(timeStamp);
 503         GetVsTimeHisto(trg.Data(),-1,-1,                                  "inV0","MultTPCvsTime")->Fill(timeStamp);
 504       }
 505
 506
 507       // has the event at least one track satisfying the required conditions? (used for rate/luminosity)
 508       if (TMath::Abs(track->Eta()) < etaCut) {
 509         // Fill histo (pt distribution with standard cuts)
 510         GetDistributionHisto(trg.Data(),kDistPt)->Fill(track->Pt());
 511         if (track->Pt() > 0.5) {
 512           atLeastPt05 = kTRUE;
 513         }
 514         if (track->Pt() > 1.0) {
 515           atLeastPt1  = kTRUE;
 516         }
 517       }
 518     }
 519
 520     // TEMPORARY LOOP: FILL DISTRIBUTION OF PT FOR 2 CLASSES OF EVENTS
 521     for (Int_t iTrack = 0; iTrack<ntracks; iTrack++) {
 522       AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTrack));
 523       // for each track
 524
 525       // track quality cuts
 526       if(!esdtrackCutsITSTPC->AcceptTrack(track)) continue;
 527
 528       // bring it to the primary vertex and compute impact parameters
 529       if(!track->RelateToVertex(vertex,fESD->GetMagneticField(),kVeryBig)) continue; // this is already done in AliReconstruction...
 530
 531       // track-to-vertex cut (see below)
 532       if(!SelectOnImpPar(track)) continue;
 533       if (TMath::Abs(track->Eta()) < etaCut){
 534         if (atLeastPt05) GetDistributionHisto(trg.Data(),kDistPt, "_atLeast05")->Fill(track->Pt());
 535         if (atLeastPt1)  GetDistributionHisto(trg.Data(),kDistPt, "_atLeast1" )->Fill(track->Pt());
 536       }
 537     }
 538     // END OF TEMPORARY LOOP
 539
 540     // Fill histos for luminosity: rate of events with at least one
 541     // track in the pseudo rapidity region |eta| < 0.8 and pt > 0.5 or
 542     // 1 GeV
 543     if (atLeastPt05) {
 544       GetVsTimeHisto(GetVsTimeHistoForLuminosityName(trg.Data(), 0.5))->Fill(timeStamp);
 545       if(fIsMC) GetEfficiencyHisto(kEffStepRec)->Fill(kEffPt05);
 546     }
 547     if (atLeastPt1)  {
 548
 549       if(GetVsTimeHisto(GetVsTimeHistoForLuminosityName(trg.Data(), 1.0))->Fill(timeStamp) < 0) {
 550         AliWarning(Form("Timestamp out of range %lld", timeStamp));
 551       };
 552       if(fIsMC) GetEfficiencyHisto(kEffStepRec)->Fill(kEffPt1);
 553     }
 554
 555
 556     // TPC TRACKS: loose selection in order to keep some BG track
 557     ntracks = fESD->GetNumberOfTracks(); // Ignore vertex selecion
 558     for(Int_t itrack = 0; itrack < ntracks; itrack++){
 559       AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(itrack));
 560       // for each track
 561       // Same as AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() , but not using DCA cut
 562
 563       AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 564
 565       esdTrackCuts->SetMinNClustersTPC(50);
 566       esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 567       esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 568
 569 //       esdTrackCuts->SetMaxDCAToVertexZ(3.2);
 570 //       esdTrackCuts->SetMaxDCAToVertexXY(2.4);
 571 //       esdTrackCuts->SetDCAToVertex2D(kTRUE);
 572
 573       if (!esdTrackCuts->AcceptTrack(track)) continue;
 574       // Fill pt and DCA distribution
 575       // pt
 576       GetDistributionHisto(trg.Data(),kDistPtLoose)->Fill(track->Pt());
 577       if(isInV0) GetDistributionHisto(trg.Data(),kDistPtLoose, "_inV0")->Fill(track->Pt());
 578       // dca (on the xy plane)
 579       Float_t xy_dca, z_dca;
 580       track->GetImpactParameters(xy_dca,z_dca);
 581       GetDistributionHisto(trg.Data(),kDistDCATPC)->Fill(xy_dca);
 582       if(isInV0) GetDistributionHisto(trg.Data(),kDistDCATPC, "_inV0")->Fill(track->Pt());
 583
 584     }
 585
 586
 587   }
 588
 589 }
 590
 591 void   AliAnalysisTaskBGvsTime::Terminate(Option_t *){
 592
 593   // normalize and rescale histos
 594   AliInfo("Normalizing multiplicity histo") ;
 595
 596   fListWrapper = dynamic_cast<AliHistoListWrapper*> (GetOutputData(1));
 597   if (!fListWrapper){
 598     AliError("Cannot get list wrapper");
 599   }
 600   fListHisto = fListWrapper->GetList();
 601   fHistoTimeStampDiffUTC = (TH1F*) fListHisto->FindObject("fHistoTimeStampDiffUTC");
 602   fHistoTimeStampVsUTC   = (TH2F*) fListHisto->FindObject("fHistoTimeStampVsUTC");
 603
 604   // Divide rate histos for the number of events
 605   TIterator * iter = fListHisto->MakeIterator();
 606   TH1 * h = 0;
 607   while ((h = (TH1*) iter->Next())) {
 608     if((TString(h->GetName())).Contains("hMultSPDvsTime_") ||
 609        (TString(h->GetName())).Contains("hMultTPCvsTime_")
 610        ){
 611
 612       AliInfo(Form("Normalizing %s",h->GetName()));
 613       //      continue;
 614       TString histoname = h->GetName();
 615       histoname.ReplaceAll("MultSPDvsTime","");
 616       histoname.ReplaceAll("MultTPCvsTime","");
 617
 618       TH1* hev = (TH1*)fListHisto->FindObject(histoname.Data());
 619       if(!hev) {
 620         AliError(Form(" -> Cannot find events histo %s",histoname.Data()));
 621         continue;
 622       }
 623       AliInfo (Form(" with histo %s",hev->GetName()));
 624       // Errors on ev num should be ignored in the division
 625       Int_t nbin  = h->GetNbinsX();
 626       for(Int_t ibin = 1; ibin <=nbin; ibin++) {
 627         if(hev->GetBinContent(ibin) != 0) {
 628           h->SetBinContent(ibin, h->GetBinContent(ibin)/ hev->GetBinContent(ibin));
 629           h->SetBinError  (ibin, h->GetBinError(ibin)  / hev->GetBinContent(ibin));
 630         }
 631       }
 632     }
 633     if((TString(h->GetName())).Contains("hRate")){
 634       h->Scale(1.,"width"); // divide for bin width to obtain a rate
 635     }
 636   }
 637
 638   // Compute Efficiency:
 639   if(fIsMC) {
 640     TH1F* hgen =(TH1F*) fListHisto->FindObject("hEffGen");
 641     TH1F* htrg =(TH1F*) fListHisto->FindObject("hEffTrig");
 642     TH1F* hrec =(TH1F*) fListHisto->FindObject("hEffRec");
 643     if (!hgen || !hrec || !htrg) {
 644       AliError("Cannot find eff histos");
 645     }
 646     else {
 647       TH1F* heff = (TH1F*) fListHisto->FindObject("hEffRatio");
 648       TH1F* hefftrg = (TH1F*) fListHisto->FindObject("hEffRatioTrg");
 649       if (heff) {
 650         AliWarning("hEffRatio already in output list?");
 651       }
 652       else {
 653         heff = (TH1F*) hgen->Clone("hEffRatio");
 654       }
 655       heff->Reset();
 656       heff->Divide(hrec,hgen,1,1,"B");
 657       fListHisto->Add(heff);
 658
 659       if (hefftrg) {
 660         AliWarning("hEffRatioTrg already in output list?");
 661       }
 662       else {
 663         hefftrg = (TH1F*) hgen->Clone("hEffRatioTrg");
 664       }
 665       hefftrg->Reset();
 666       hefftrg->Divide(htrg,hgen,1,1,"B");
 667       fListHisto->Add(hefftrg);
 668
 669     }
 670   }
 671
 672
 673   AliInfo(Form("Time interval: %lld -- %lld",fFirstTimeStamp,fLastTimeStamp));
 674
 675   AliInfo("Saving physics selection histos");
 676   fPhysicsSelection = dynamic_cast<AliPhysicsSelection*> (GetOutputData(2));
 677
 678   TFile* fout = new TFile("event_stat.root", "RECREATE");
 679
 680   if (fPhysicsSelection)
 681       {
 682         fPhysicsSelection->Print();
 683         fPhysicsSelection->SaveHistograms("physics_selection");
 684       }
 685
 686   fout->Write();
 687   fout->Close();
 688
 689 }
 690
 691 TH1F * AliAnalysisTaskBGvsTime::BookVsTimeHisto(const char * name, const char * title){
 692
 693   // Book istograms vs time
 694
 695   AliInfo(Form("Booking histo %s",name));
 696
 697   Bool_t oldStatus = TH1::AddDirectoryStatus();
 698   TH1::AddDirectory(kFALSE);
 699
 700   //  static Int_t nBins =  fEndTime - fStartTime + 1;
 701   // Compute bin width and range.
 702   // if time interval is not a multiple of bin width drop the last bin
 703
 704   Float_t lenght =  fEndTime - fStartTime;
 705   Int_t nBins = TMath::FloorNint(lenght/fBinWidth)+1;
 706   Float_t * bins = new Float_t[nBins+1];
 707   for(Int_t ibin = 0; ibin <= nBins; ibin++){
 708     Float_t edge = ibin*fBinWidth - 0.5;
 709     bins[ibin] = edge < (lenght - 0.5) ? edge  : (lenght-0.5);
 710   }
 711
 712   //  Int_t end   = nBins*fBinWidth;
 713
 714   //  TH1F * h = new TH1F(name,title, nBins,  - 0.5, fEndTime - fStartTime + 0.5);
 715   TH1F * h = new TH1F(name,title, nBins, bins);
 716   fListHisto->Add(h);
 717   h->SetXTitle("orbit");
 718   h->Sumw2();
 719
 720   TH1::AddDirectory(oldStatus);
 721
 722
 723   delete [] bins;
 724   return h;
 725 }
 726
 727 TH1F * AliAnalysisTaskBGvsTime::GetVsTimeHisto(const char * name) {
 728
 729   // Returns vs time histo. If not existing, creates it.
 730
 731   TH1F * h = (TH1F*) fListHisto->FindObject(name);
 732   if(!h) h = BookVsTimeHisto(name,name);
 733   return h;
 734
 735 }
 736
 737 TH1F * AliAnalysisTaskBGvsTime::GetEfficiencyHisto(Int_t step) {
 738
 739   // Return efficiency histo. If not existing, creates it.
 740   // 1 bin per category
 741
 742   // Probability of reconstructing one event with at least one
 743   // particle in a given category
 744   // This is only created for MC.
 745
 746   // If the first argument is true, returns the histo at generation
 747   // level, otherwise returns the histo at rec level. Fill this histo
 748   // with elements of the enum kEffPt1, kEffPt05, ...
 749
 750   const char * name = 0;
 751   switch(step) {
 752   case kEffStepGen:
 753     name =  "hEffGen";
 754     break;
 755   case kEffStepRec:
 756     name =  "hEffRec";
 757     break;
 758   case kEffStepTrig:
 759     name =  "hEffTrig";
 760     break;
 761   }
 762
 763
 764   TH1F * h = (TH1F*) fListHisto->FindObject(name);
 765   if(!h) {
 766     Bool_t oldStatus = TH1::AddDirectoryStatus();
 767     TH1::AddDirectory(kFALSE);
 768     h = new TH1F (name,name, kNEff, -0.5, kNEff-0.5);
 769     h->GetXaxis()->SetBinLabel(h->FindBin(kEffPt1) , "pt > 1.0 GeV");
 770     h->GetXaxis()->SetBinLabel(h->FindBin(kEffPt05), "pt > 0.5 GeV");
 771     h->Sumw2();
 772     fListHisto->Add(h);
 773     TH1::AddDirectory(oldStatus);
 774   }
 775
 776   return h;
 777
 778 }
 779
 780 const char * AliAnalysisTaskBGvsTime::GetVsTimeHistoForLuminosityName(const char * triggerClass, Float_t ptmin) {
 781
 782   // Compose name of rate histos
 783
 784   static TString name;
 785   name = "hRate_";
 786   name += triggerClass;
 787
 788   name += Form ("_etamin_0.8_ptmin_%2.2f",ptmin);
 789
 790   return name.Data();
 791 }
 792
 793
 794 const char * AliAnalysisTaskBGvsTime::GetVsTimeHistoName(const char * triggerClass, Int_t nclusters, Int_t bx, const char * v0flag, const char *prefix){
 795
 796   // compose name of vs time histos w/ different cuts
 797
 798   static TString name;
 799   name = "h";
 800   name = name +prefix+"_";
 801   name += triggerClass;
 802
 803   if (nclusters >= 0) {
 804     if(fMultBins){
 805       // Add multiplicity label
 806       Int_t selected_bin = -1;
 807       for(Int_t ibin = 0; ibin < (fNMultBins-1); ibin++){
 808         if (nclusters < fMultBins[ibin+1]) {
 809           selected_bin = ibin;
 810           break;
 811         }
 812       }
 813       if(selected_bin < 0) {
 814         AliError(Form("Cannot get mult bin - %d", nclusters));
 815         selected_bin = (fNMultBins-1);
 816       }
 817       name += Form("_SPDCls_%d-%d",fMultBins[selected_bin],fMultBins[selected_bin+1]);
 818       //    AliInfo(Form("Mult bin: %d, %d, %s", nclusters, selected_bin, name.Data()));
 819     }
 820   } else {
 821     name += "_SPDCls_ALL";
 822   }
 823   if(bx < 0) name = name +"_BXALL_"+v0flag;
 824   else       name = name +"_BX"+long(bx)+"_"+v0flag;
 825   return name.Data();
 826 }
 827 const char * AliAnalysisTaskBGvsTime::GetVsTimeHistoNameAll(const char * triggerClass) {
 828
 829   // Compose default name of vstime histos in a given trigger class
 830
 831   static TString name;
 832   name = "h_";
 833   name += triggerClass;
 834   return name.Data();
 835
 836 }
 837
 838 TH1F * AliAnalysisTaskBGvsTime::GetDistributionHisto(const char * triggerClass, Int_t dist, const char * suffix) {
 839
 840   // Returns distributions histos. If not existing, creates it.
 841
 842   // Possible distributions:
 843   // - SPD cluster
 844   // - TPC tracks multiplicity
 845   // - TPC tracks pt
 846   // - Vz
 847   // - TPC tracks DCA
 848   // - cluster per ITS layer
 849
 850   TString name  ;
 851   TString title ;
 852   TString xtitle;
 853   Int_t   nbin =0;
 854   Float_t   min  =0;
 855   Float_t   max  =0;
 856
 857
 858   if (dist == kDistSPDMult) {
 859     name  = "hMultSPD_";
 860     title = "SPD Cluster Multiplicity (";
 861     xtitle = "SPD Clusters";
 862     nbin  = 50;
 863     min   = 0;
 864     max   = 100;
 865   } else if (dist == kDistTPCMult) {
 866     name  = "hMultTPC_";
 867     title = "TPC Tracks Multiplicity (";
 868     xtitle = "TPC Tracks";
 869     nbin  = 25;
 870     min   = 0;
 871     max   = 50;
 872   } else if (dist == kDistPtLoose) {
 873     name  = "hPtLoose_";
 874     title = "p_{T} distribution - TPC, loose cuts (";
 875     xtitle = "p_{T} (GeV)";
 876     nbin  = 50;
 877     min   = 0;
 878     max   = 10;
 879   } else if (dist == kDistPt) {
 880     name  = "hPt_";
 881     title = "p_{T} distribution - TPC, standard cuts, |#eta|<0.8 (";
 882     xtitle = "p_{T} (GeV)";
 883     nbin  = 100;
 884     min   = 0;
 885     max   = 10;
 886   } else if (dist == kDistVertex || dist == kDistVertexZ || dist == kDistVertex3D ) {
 887     nbin  = 120;
 888     min   = -30;
 889     max   = 30;
 890     name  = "hVz_";
 891     title = "V_{z} distribution";
 892     if      (dist == kDistVertexZ)  {
 893       title += " - Vertexer Z  (";
 894       name+= "Z_";
 895     }
 896     else if (dist == kDistVertex3D) {
 897       title += " - Vertexer 3D (";
 898       name+= "3D_";
 899     }
 900     else    title += " (";
 901     xtitle = "V_{z} (cm)";
 902   } else if (dist == kDistDCATPC) {
 903     nbin  = 50;
 904     min   = -0.5;
 905     max   = +0.5;
 906     name  = "hDCA_";
 907     title = "TPC DCA distribution - loose cuts (";
 908     xtitle = "DCA";
 909   } else if (dist == kDistClsITSLayer) {
 910     nbin  = 6;
 911     min   = -0.5;
 912     max   = 5.5;
 913     name  = "hClsITS_";
 914     title = "Cluster per ITS layer (";
 915     xtitle = "ITS Layer";
 916   }else {
 917     AliError(Form("Distribution type not supported: %d",dist));
 918     return 0;
 919   }
 920
 921
 922   name += triggerClass;
 923   if (suffix) name += suffix;
 924   title = title+name+")";
 925
 926   TH1F* h = (TH1F*) fListHisto->FindObject(name.Data());
 927
 928   if(!h) h = BookDistributionHisto(name.Data(), title.Data(), xtitle.Data(), nbin, min, max);
 929   return h;
 930
 931 }
 932
 933 TH1F * AliAnalysisTaskBGvsTime::BookDistributionHisto(const char * name, const char * title, const char * xtitle, Int_t nbin, Float_t min, Float_t max) {
 934
 935   // Book distributions histos
 936
 937   AliInfo(Form("Booking histo %s",name));
 938
 939   Bool_t oldStatus = TH1::AddDirectoryStatus();
 940   TH1::AddDirectory(kFALSE);
 941
 942   TH1F * h = new TH1F (name,title,nbin,min,max);
 943   h->Sumw2();
 944   h->SetXTitle(xtitle);
 945   //  h->SetYTitle("N");
 946   fListHisto->Add(h);
 947   TH1::AddDirectory(oldStatus);
 948
 949   return h;
 950
 951 }
 952
 953 TH1F * AliAnalysisTaskBGvsTime::GetDeadTimeHisto(const char * triggerClass) {
 954   // returns histo of dead time vs timestamp for a given trigger
 955   // class. If the histo does not exist, it books it.
 956   TString name  = "hDeadTime_"  ;
 957   TString title = "Dead Time vs TimeStamp (";
 958
 959   name += triggerClass;
 960   title = title+triggerClass+")";
 961
 962
 963   TH1F* h = (TH1F*) fListHisto->FindObject(name.Data());
 964
 965   if(!h) {
 966     h = BookVsTimeHisto(name.Data(), title.Data());
 967     h->SetYTitle("deadtime");
 968   }
 969   return h;
 970
 971
 972 }
 973
 974
 975 // TH2F * AliAnalysisTaskBGvsTime::BookDeadTimeHisto(const char * name, const char * title) {
 976 //   // Book dead time vs time stamp histos
 977
 978 //   AliInfo(Form("Booking histo %s",name));
 979
 980 //   Bool_t oldStatus = TH1::AddDirectoryStatus();
 981 //   TH1::AddDirectory(kFALSE);
 982
 983 //   Float_t lenght =  fEndTime - fStartTime;
 984 //   Int_t nBins = TMath::FloorNint(lenght/fBinWidth)+1;
 985 //   Double_t * bins = new Double_t[nBins+1];
 986 //   for(Int_t ibin = 0; ibin <= nBins; ibin++){
 987 //     Float_t edge = ibin*fBinWidth - 0.5;
 988 //     bins[ibin] = edge < (lenght - 0.5) ? edge  : (lenght-0.5);
 989 //   }
 990
 991 //   TH2F * h = new TH2F (name,title,nBins,bins,200,0.,2.);
 992 //   h->Sumw2();
 993 //   h->SetXTitle("Time (s)");
 994 //   h->SetYTitle("Dead Time");
 995 //   fListHisto->Add(h);
 996 //   TH1::AddDirectory(oldStatus);
 997
 998 //   delete bins;
 999
1000 //   return h;
1001
1002 // }
1003
1004
1005
1006 Bool_t AliAnalysisTaskBGvsTime::IsEventInBinZero() {
1007
1008   // Returns true if an event is to be assigned to the zero bin
1009
1010   Bool_t isZeroBin = kTRUE;
1011   const AliESDEvent* esd=fESD;
1012   const AliMultiplicity* mult = esd->GetMultiplicity();
1013   if (!mult){
1014     Printf("AliAnalysisTaskBGvsTime::IsBinZero: Can't get mult object");
1015     return kFALSE;
1016   }
1017   Int_t ntracklet = mult->GetNumberOfTracklets();
1018   const AliESDVertex * vtxESD = esd->GetPrimaryVertexSPD();
1019   if(vtxESD) {
1020     // If there is a vertex from vertexer z with delta phi > 0.02 we
1021     // don't consider it rec (we keep the event in bin0). If quality
1022     // is good eneough we check the number of tracklets
1023     // if the vertex is more than 15 cm away, this is autamatically bin0
1024     if( TMath::Abs(vtxESD->GetZ()) <= 15 ) {
1025       if (vtxESD->IsFromVertexerZ()) {
1026         if (vtxESD->GetDispersion()<=0.02 ) {
1027           if(ntracklet>0) isZeroBin = kFALSE;
1028         }
1029       } else if(ntracklet>0) isZeroBin = kFALSE; // if the event is not from Vz we chek the n of tracklets
1030     }
1031   }
1032   return isZeroBin;
1033
1034 }
1035
1036 Bool_t AliAnalysisTaskBGvsTime::SelectOnImpPar(AliESDtrack* t) {
1037   // from andrea dainese
1038   // cut on transverse impact parameter
1039   Float_t d0z0[2],covd0z0[3];
1040   t->GetImpactParameters(d0z0,covd0z0);
1041   Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
1042   Float_t d0max = 7.*sigma;
1043   if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
1044   return kFALSE;
1045 }
1046
1047
1048
1049