PWGUD/selectors/multiplicity/AliMultiplicityTask.cxx

   1 /* $Id$ */
   2
   3 #include "AliMultiplicityTask.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TVector3.h>
   9 #include <TChain.h>
  10 #include <TFile.h>
  11 #include <TH1D.h>
  12 #include <TH2F.h>
  13 #include <TH3F.h>
  14 #include <TParticle.h>
  15 #include <TRandom.h>
  16 #include <TNtuple.h>
  17 #include <TObjString.h>
  18 #include <TF1.h>
  19
  20 #include <AliLog.h>
  21 #include <AliESDVertex.h>
  22 #include <AliESDEvent.h>
  23 #include <AliStack.h>
  24 #include <AliHeader.h>
  25 #include <AliGenEventHeader.h>
  26 #include <AliMultiplicity.h>
  27 #include <AliAnalysisManager.h>
  28 #include <AliMCEventHandler.h>
  29 #include <AliMCEvent.h>
  30 #include <AliESDInputHandler.h>
  31
  32 #include "AliESDtrackCuts.h"
  33 #include "AliPWG0Helper.h"
  34 #include "AliMultiplicityCorrection.h"
  35 #include "AliCorrection.h"
  36 #include "AliCorrectionMatrix3D.h"
  37 #include "AliPhysicsSelection.h"
  38 #include "AliTriggerAnalysis.h"
  39 #include "AliVEvent.h"
  40
  41 ClassImp(AliMultiplicityTask)
  42
  43 AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
  44   AliAnalysisTask("AliMultiplicityTask", ""),
  45   fESD(0),
  46   fOption(opt),
  47   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn)),
  48   fTrigger(AliTriggerAnalysis::kMB1),
  49   fDeltaPhiCut(-1),
  50   fDiffTreatment(AliPWG0Helper::kMCFlags),
  51   fReadMC(kFALSE),
  52   fUseMCVertex(kFALSE),
  53   fMultiplicity(0),
  54   fEsdTrackCuts(0),
  55   fSystSkipParticles(kFALSE),
  56   fSelectProcessType(0),
  57   fParticleSpecies(0),
  58   fdNdpT(0),
  59   fPtSpectrum(0),
  60   fTemp1(0),
  61   fTemp2(0),
  62   fVertex(0),
  63   fEtaPhi(0),
  64   fOutput(0)
  65 {
  66   //
  67   // Constructor. Initialization of pointers
  68   //
  69
  70   for (Int_t i = 0; i<8; i++)
  71     fParticleCorrection[i] = 0;
  72
  73   for (Int_t i=0; i<3; i++)
  74     fEta[i] = 0;
  75
  76   // Define input and output slots here
  77   DefineInput(0, TChain::Class());
  78   DefineOutput(0, TList::Class());
  79
  80   if (fOption.Contains("only-process-type-nd"))
  81     fSelectProcessType = 1;
  82
  83   if (fOption.Contains("only-process-type-sd"))
  84     fSelectProcessType = 2;
  85
  86   if (fOption.Contains("only-process-type-dd"))
  87     fSelectProcessType = 3;
  88
  89   if (fSelectProcessType != 0)
  90     AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
  91 }
  92
  93 AliMultiplicityTask::~AliMultiplicityTask()
  94 {
  95   //
  96   // Destructor
  97   //
  98
  99   // histograms are in the output list and deleted when the output
 100   // list is deleted by the TSelector dtor
 101
 102   if (fOutput&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 103     delete fOutput;
 104     fOutput = 0;
 105   }
 106 }
 107
 108 //________________________________________________________________________
 109 void AliMultiplicityTask::ConnectInputData(Option_t *)
 110 {
 111   // Connect ESD
 112   // Called once
 113
 114   Printf("AliMultiplicityTask::ConnectInputData called");
 115
 116   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 117   if (!tree) {
 118     Printf("ERROR: Could not read tree from input slot 0");
 119   } else {
 120     // Disable all branches and enable only the needed ones
 121     /*
 122     tree->SetBranchStatus("*", 0);
 123
 124     tree->SetBranchStatus("AliESDHeader*", 1);
 125     tree->SetBranchStatus("*Vertex*", 1);
 126
 127     if (fAnalysisMode & AliPWG0Helper::kSPD) {
 128       tree->SetBranchStatus("AliMultiplicity*", 1);
 129     }
 130
 131     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
 132       //AliESDtrackCuts::EnableNeededBranches(tree);
 133       tree->SetBranchStatus("Tracks*", 1);
 134     }
 135     */
 136
 137     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 138
 139     if (!esdH) {
 140       Printf("ERROR: Could not get ESDInputHandler");
 141     } else
 142       fESD = esdH->GetEvent();
 143   }
 144
 145   // disable info messages of AliMCEvent (per event)
 146   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 147 }
 148
 149 void AliMultiplicityTask::CreateOutputObjects()
 150 {
 151   // create result objects and add to output list
 152
 153   fOutput = new TList;
 154   fOutput->SetOwner();
 155
 156   fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 157   fOutput->Add(fMultiplicity);
 158
 159   fdNdpT = new TH1F("fdNdpT", "fdNdpT", 1000, 0, 10);
 160   fdNdpT->Sumw2();
 161   fOutput->Add(fdNdpT);
 162
 163   if (fOption.Contains("particle-efficiency"))
 164     for (Int_t i = 0; i<8; i++)
 165     {
 166       fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
 167       fOutput->Add(fParticleCorrection[i]);
 168     }
 169
 170   if (fOption.Contains("pt-spectrum-hist"))
 171   {
 172     TFile* file = TFile::Open("ptspectrum_fit.root");
 173     if (file)
 174     {
 175       TString subStr(fOption(fOption.Index("pt-spectrum")+17, 3));
 176       TString histName(Form("ptspectrum_%s", subStr.Data()));
 177       AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
 178       fPtSpectrum = (TH1D*) file->Get(histName);
 179       if (!fPtSpectrum)
 180         AliError("Histogram not found");
 181     }
 182     else
 183       AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
 184   }
 185
 186   if (fOption.Contains("pt-spectrum-func"))
 187   {
 188     if (fPtSpectrum)
 189     {
 190       Printf("Using function for systematic p_t study");
 191     }
 192     else
 193     {
 194       Printf("ERROR: Could not find function for systematic p_t study");
 195       fPtSpectrum = new TH1D("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
 196       fPtSpectrum->SetBinContent(1, 1);
 197     }
 198   }
 199
 200   if (fPtSpectrum)
 201     Printf("WARNING: Systematic study enabled. Pt spectrum will be modified");
 202
 203   if (fOption.Contains("particle-species")) {
 204     fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "vtx:Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec:nolabel:doublePrim:doubleCount");
 205     fOutput->Add(fParticleSpecies);
 206   }
 207
 208   fTemp1 = new TH2F("fTemp1", "fTemp1", 100, -0.5, 99.5, 100, -0.5, 99.5);
 209   fOutput->Add(fTemp1);
 210
 211   for (Int_t i=0; i<3; i++)
 212   {
 213     fEta[i] = new TH1F(Form("fEta_%d", i), ";#eta", 100, -2, 2);
 214     fOutput->Add(fEta[i]);
 215   }
 216
 217   fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
 218   fOutput->Add(fVertex);
 219
 220   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
 221   fOutput->Add(fEtaPhi);
 222
 223   // TODO set seed for random generator
 224 }
 225
 226 void AliMultiplicityTask::Exec(Option_t*)
 227 {
 228   // process the event
 229
 230   // Check prerequisites
 231   if (!fESD)
 232   {
 233     AliDebug(AliLog::kError, "ESD not available");
 234     return;
 235   }
 236
 237   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 238   if (!inputHandler)
 239   {
 240     Printf("ERROR: Could not receive input handler");
 241     return;
 242   }
 243
 244   Bool_t eventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
 245
 246   static AliTriggerAnalysis* triggerAnalysis = 0;
 247   if (!triggerAnalysis)
 248   {
 249     AliPhysicsSelection* physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 250     triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 251   }
 252   if (eventTriggered)
 253     eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
 254
 255   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 256   if (vtxESD && !AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
 257     vtxESD = 0;
 258
 259   // remove vertices outside +- 15 cm
 260   if (vtxESD && TMath::Abs(vtxESD->GetZ()) > 15)
 261     vtxESD = 0;
 262
 263   Bool_t eventVertex = (vtxESD != 0);
 264
 265   Double_t vtx[3];
 266   if (vtxESD)
 267   {
 268     vtxESD->GetXYZ(vtx);
 269     fVertex->Fill(vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ());
 270   }
 271
 272   // post the data already here
 273   PostData(0, fOutput);
 274
 275   //const Float_t kPtCut = 0.3;
 276
 277   // create list of (label, eta) tuples
 278   Int_t inputCount = 0;
 279   Int_t* labelArr = 0;
 280   Float_t* etaArr = 0;
 281   Int_t nGoodTracks = 0;  // number of total good tracks is needed both in SPD and TPC loops if the mode is kTPCSPD
 282   Int_t nESDTracks  = 0; // Total number of ESD tracks (including those not passing the cuts);
 283   const int kRejBit = BIT(15); // set this bit in ESD tracks if it is rejected by a cut
 284
 285   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS || fAnalysisMode & AliPWG0Helper::kTPCSPD)
 286   {
 287     if (!fEsdTrackCuts)
 288     {
 289       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 290       return;
 291     }
 292
 293     if (vtxESD)
 294     {
 295       // get multiplicity from ESD tracks
 296       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
 297       nGoodTracks = list->GetEntries();
 298
 299       Int_t arraySize = nGoodTracks; // if we're also using tracklets, we need to increase this
 300       if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
 301         // if I have to also count clusters later on, I have to make sure the arrays are big enough
 302         // Moreover, we need to keep track of rejected tracks to avoid double-counting
 303         Printf( "Processing tracks");
 304         const AliMultiplicity* mult = fESD->GetMultiplicity();
 305         if (!mult)
 306           {
 307             AliDebug(AliLog::kError, "AliMultiplicity not available");
 308             return;
 309           }
 310         arraySize += mult->GetNumberOfTracklets();
 311         // allocate and fill array for rejected tracks
 312         nESDTracks = fESD->GetNumberOfTracks();
 313         for(Int_t itracks = 0; itracks < nESDTracks; itracks++){
 314           AliESDtrack* track = fESD->GetTrack(itracks);
 315           if( itracks!=track->GetID() ){
 316             AliFatal("Track ID not matching track index!");
 317           }
 318           if (!fEsdTrackCuts->AcceptTrack(track)) {
 319             track->SetBit(kRejBit);
 320           }
 321         }
 322       }
 323
 324       labelArr = new Int_t[arraySize];
 325       etaArr = new Float_t[arraySize];
 326
 327       // loop over esd tracks
 328       for (Int_t i=0; i<nGoodTracks; i++)
 329       {
 330         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 331         if (!esdTrack)
 332         {
 333           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 334           continue;
 335         }
 336
 337         if (esdTrack->Pt() < 0.15){
 338           esdTrack->SetBit(kRejBit);
 339           continue;
 340         }
 341
 342         if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
 343           // Cuts by ruben to flag secondaries
 344           if (esdTrack->IsOn(AliESDtrack::kMultSec) ) continue;
 345           if (esdTrack->IsOn(AliESDtrack::kMultInV0)) continue;
 346         }
 347
 348         Float_t d0z0[2],covd0z0[3];
 349         esdTrack->GetImpactParameters(d0z0,covd0z0);
 350         Float_t sigma= 0.0050+0.0060/TMath::Power(esdTrack->Pt(),0.9);
 351         Float_t d0max = 7.*sigma;
 352         if (TMath::Abs(d0z0[0]) > d0max) {
 353           //      rejLabelArr[irejCount++] = esdTrack->GetID(); // We
 354           //      don't count the tracklet if the track is a
 355           //      secondary, so this must be commented out
 356           continue;
 357         }
 358
 359         if (vtxESD && TMath::Abs(vtx[2]) < 10)
 360           fEtaPhi->Fill(esdTrack->Eta(), esdTrack->Phi());
 361
 362         etaArr[inputCount] = esdTrack->Eta();
 363         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 364         ++inputCount;
 365       }
 366
 367       delete list;
 368     }
 369   }
 370   // SPD tracklets are used either in SPD standalone mode or in TPC+SPD mode (tracks + tracklets not using any cluster already used in tracks)
 371   if (fAnalysisMode & AliPWG0Helper::kSPD || fAnalysisMode & AliPWG0Helper::kTPCSPD)
 372   {
 373
 374     if (vtxESD)
 375     {
 376       AliDebug(AliLog::kInfo, "Processing tracklets");
 377       // get tracklets
 378       const AliMultiplicity* mult = fESD->GetMultiplicity();
 379       if (!mult)
 380       {
 381         AliDebug(AliLog::kError, "AliMultiplicity not available");
 382         if (labelArr) delete[] labelArr;
 383         if (etaArr) delete[] etaArr;
 384         return;
 385       }
 386
 387       Bool_t foundInEta10 = kFALSE;
 388       if ((fAnalysisMode & AliPWG0Helper::kSPD) && !(fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS || fAnalysisMode & AliPWG0Helper::kTPCSPD)) {
 389         // if we are counting both tracks and tracklets, these arrays were already initialized above
 390         AliDebug(AliLog::kInfo, "Booking arrays");
 391         if (!labelArr) labelArr = new Int_t[mult->GetNumberOfTracklets()];
 392         if (!etaArr) etaArr = new Float_t[mult->GetNumberOfTracklets()];
 393       }
 394       if (inputCount) foundInEta10 = kTRUE; // by definition, if we found a track.
 395
 396       // get multiplicity from ITS tracklets
 397       for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 398       {
 399         //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 400
 401         Float_t deltaPhi = mult->GetDeltaPhi(i);
 402
 403         if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
 404           continue;
 405
 406         if (fSystSkipParticles && gRandom->Uniform() < (0.0153))
 407           {
 408           Printf("Skipped tracklet!");
 409           continue;
 410           }
 411
 412         // if counting tracks+tracklets, check if clusters where already used in tracks
 413         if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
 414             Int_t id1,id2;
 415             mult->GetTrackletTrackIDs(i,0,id1,id2);
 416             if ( (id1>=0&& !fESD->GetTrack(id1)->TestBit(kRejBit) ) ||
 417                  (id2>=0&& !fESD->GetTrack(id2)->TestBit(kRejBit) )
 418                  ) {
 419               printf ("Skipping tracklet: already used in track");
 420               continue;
 421             }
 422         // Ruben also had this, but we're not using pure ITSSA tracks here:
 423         // if (fSPDMult->FreeClustersTracklet(itr,1)) trITSSApure++; // not used in ITS_SA_Pure track
 424         }
 425
 426
 427         etaArr[inputCount] = mult->GetEta(i);
 428         if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
 429         {
 430           labelArr[inputCount] = mult->GetLabel(i, 0);
 431         }
 432         else
 433           labelArr[inputCount] = -1;
 434
 435         for (Int_t j=0; j<3; j++)
 436         {
 437           if (vtx[2] > fMultiplicity->GetVertexBegin(j) && vtx[2] < fMultiplicity->GetVertexEnd(j))
 438             fEta[j]->Fill(etaArr[inputCount]);
 439         }
 440
 441         if (vtxESD && TMath::Abs(vtx[2]) < 10)
 442           fEtaPhi->Fill(etaArr[inputCount], mult->GetPhi(i));
 443
 444           // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
 445         if (TMath::Abs(etaArr[inputCount]) < 1)
 446           foundInEta10 = kTRUE;
 447
 448         ++inputCount;
 449       }
 450
 451       if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
 452         eventTriggered = kFALSE;
 453     }
 454   }
 455
 456   // eta range for nMCTracksSpecies, nESDTracksSpecies
 457   Float_t etaRange = 1.0;
 458 //   switch (fAnalysisMode) {
 459 //     case AliPWG0Helper::kInvalid: break;
 460 //     case AliPWG0Helper::kTPC:
 461 //     case AliPWG0Helper::kTPCITS:
 462 //      etaRange = 1.0; break;
 463 //     case AliPWG0Helper::kSPD: etaRange = 1.0; break;
 464 //     default: break;
 465 //   }
 466
 467   if (!fReadMC) // Processing of ESD information
 468   {
 469     Int_t nESDTracks05 = 0;
 470     Int_t nESDTracks10 = 0;
 471     Int_t nESDTracks14 = 0;
 472
 473     for (Int_t i=0; i<inputCount; ++i)
 474     {
 475       Float_t eta = etaArr[i];
 476
 477       if (TMath::Abs(eta) < 0.5)
 478         nESDTracks05++;
 479
 480       if (TMath::Abs(eta) < 1.0)
 481         nESDTracks10++;
 482
 483       if (TMath::Abs(eta) < 1.3)
 484         nESDTracks14++;
 485     }
 486
 487     //if (nESDTracks05 >= 20 || nESDTracks10 >= 30 || nESDTracks14 >= 32)
 488     //  Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d; Tracks: %d %d %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber(), nESDTracks05, nESDTracks10, nESDTracks14);
 489
 490     if (eventTriggered)
 491       fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
 492
 493     if (eventTriggered && eventVertex)
 494       fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 495   }
 496   else if (fReadMC)   // Processing of MC information
 497   {
 498     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 499     if (!eventHandler) {
 500       Printf("ERROR: Could not retrieve MC event handler");
 501       if (labelArr) delete[] labelArr;
 502       if (etaArr) delete[] etaArr;
 503       return;
 504     }
 505
 506     AliMCEvent* mcEvent = eventHandler->MCEvent();
 507     if (!mcEvent) {
 508       Printf("ERROR: Could not retrieve MC event");
 509       if (labelArr) delete[] labelArr;
 510       if (etaArr) delete[] etaArr;
 511       return;
 512     }
 513
 514     AliStack* stack = mcEvent->Stack();
 515     if (!stack)
 516     {
 517       AliDebug(AliLog::kError, "Stack not available");
 518       if (labelArr) delete[] labelArr;
 519       if (etaArr) delete[] etaArr;
 520       return;
 521     }
 522
 523     AliHeader* header = mcEvent->Header();
 524     if (!header)
 525     {
 526       AliDebug(AliLog::kError, "Header not available");
 527       if (labelArr) delete[] labelArr;
 528       if (etaArr) delete[] etaArr;
 529       return;
 530     }
 531
 532     if (fUseMCVertex)
 533     {
 534       Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 535       // get the MC vertex
 536       AliGenEventHeader* genHeader = header->GenEventHeader();
 537       TArrayF vtxMC(3);
 538       genHeader->PrimaryVertex(vtxMC);
 539
 540       vtx[2] = vtxMC[2];
 541     }
 542
 543     // get process information
 544     AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
 545
 546     Bool_t processEvent = kTRUE;
 547     if (fSelectProcessType > 0)
 548     {
 549       processEvent = kFALSE;
 550
 551       // non diffractive
 552       if (fSelectProcessType == 1 && processType == AliPWG0Helper::kND)
 553         processEvent = kTRUE;
 554
 555       // single diffractive
 556       if (fSelectProcessType == 2 && processType == AliPWG0Helper::kSD)
 557         processEvent = kTRUE;
 558
 559       // double diffractive
 560       if (fSelectProcessType == 3 && processType == AliPWG0Helper::kDD)
 561         processEvent = kTRUE;
 562
 563       if (!processEvent)
 564         Printf("Skipping this event, because it is not of the requested process type (%d)", (Int_t) processType);
 565     }
 566
 567     if (processEvent)
 568     {
 569       // get the MC vertex
 570       AliGenEventHeader* genHeader = header->GenEventHeader();
 571       TArrayF vtxMC(3);
 572       genHeader->PrimaryVertex(vtxMC);
 573
 574       // get number of tracks from MC
 575       // loop over mc particles
 576       Int_t nPrim  = stack->GetNprimary();
 577       Int_t nMCPart = stack->GetNtrack();
 578
 579       // tracks in different eta ranges
 580       Int_t nMCTracks05 = 0;
 581       Int_t nMCTracks10 = 0;
 582       Int_t nMCTracks14 = 0;
 583       Int_t nMCTracksAll = 0;
 584
 585       // tracks per particle species, in |eta| < 2 (systematic study)
 586       Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
 587       for (Int_t i = 0; i<4; ++i)
 588         nMCTracksSpecies[i] = 0;
 589
 590       for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 591       {
 592         AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 593
 594         TParticle* particle = stack->Particle(iMc);
 595
 596         if (!particle)
 597         {
 598           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 599           continue;
 600         }
 601
 602         Bool_t debug = kFALSE;
 603         if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
 604         {
 605           //printf("%d) DROPPED ", iMC);
 606           //particle->Print();
 607           continue;
 608         }
 609
 610         //printf("%d) OK      ", iMC);
 611         //particle->Print();
 612
 613         //if (particle->Pt() < kPtCut)
 614         //  continue;
 615
 616         Int_t particleWeight = 1;
 617
 618         Float_t pt = particle->Pt();
 619
 620         // in case of systematic study, weight according to the change of the pt spectrum
 621         // it cannot be just multiplied because we cannot count "half of a particle"
 622         // instead a random generator decides if the particle is counted twice (if value > 1)
 623         // or not (if value < 0)
 624         if (fPtSpectrum)
 625         {
 626           Int_t bin = fPtSpectrum->FindBin(pt);
 627           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 628           {
 629             Float_t factor = fPtSpectrum->GetBinContent(bin);
 630             if (factor > 0)
 631             {
 632               Float_t random = gRandom->Uniform();
 633               if (factor > 1 && random < factor - 1)
 634               {
 635                 particleWeight = 2;
 636               }
 637               else if (factor < 1 && random < 1 - factor)
 638                 particleWeight = 0;
 639             }
 640           }
 641         }
 642
 643         //Printf("MC weight is: %d", particleWeight);
 644
 645         if (TMath::Abs(particle->Eta()) < 0.5)
 646           nMCTracks05 += particleWeight;
 647
 648         if (TMath::Abs(particle->Eta()) < 1.0)
 649           nMCTracks10 += particleWeight;
 650
 651         if (TMath::Abs(particle->Eta()) < 1.3)
 652           nMCTracks14 += particleWeight;
 653
 654         nMCTracksAll += particleWeight;
 655
 656         if (particle->Pt() > 0 && TMath::Abs(particle->Eta()) < 1.0)
 657           fdNdpT->Fill(particle->Pt(), 1.0 / particle->Pt());
 658
 659         if (fParticleCorrection[0] || fParticleSpecies)
 660         {
 661           Int_t id = -1;
 662           switch (TMath::Abs(particle->GetPdgCode()))
 663           {
 664             case 211: id = 0; break;
 665             case 321: id = 1; break;
 666             case 2212: id = 2; break;
 667             default: id = 3; break;
 668           }
 669
 670           if (TMath::Abs(particle->Eta()) < etaRange)
 671             nMCTracksSpecies[id]++;
 672
 673           if (fParticleCorrection[id])
 674             fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
 675         }
 676       } // end of mc particle
 677
 678       fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
 679
 680       // ESD processing
 681       Int_t nESDTracks05 = 0;
 682       Int_t nESDTracks10 = 0;
 683       Int_t nESDTracks14 = 0;
 684
 685       // tracks per particle species, in |eta| < 2 (systematic study)
 686       Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
 687       for (Int_t i = 0; i<7; ++i)
 688         nESDTracksSpecies[i] = 0;
 689
 690       Bool_t* foundPrimaries = new Bool_t[nMCPart];   // to prevent double counting
 691       for (Int_t i=0; i<nPrim; i++)
 692         foundPrimaries[i] = kFALSE;
 693
 694       Bool_t* foundPrimaries2 = new Bool_t[nMCPart];   // to prevent double counting
 695       for (Int_t i=0; i<nPrim; i++)
 696         foundPrimaries2[i] = kFALSE;
 697
 698       Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
 699       for (Int_t i=0; i<nMCPart; i++)
 700         foundTracks[i] = kFALSE;
 701
 702       for (Int_t i=0; i<inputCount; ++i)
 703       {
 704         Float_t eta = etaArr[i];
 705         Int_t label = labelArr[i];
 706
 707         Int_t particleWeight = 1;
 708
 709         // in case of systematic study, weight according to the change of the pt spectrum
 710         if (fPtSpectrum)
 711         {
 712           TParticle* mother = 0;
 713
 714           // preserve label for later
 715           Int_t labelCopy = label;
 716           if (labelCopy >= 0)
 717             labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
 718           if (labelCopy >= 0)
 719             mother = stack->Particle(labelCopy);
 720
 721           // in case of pt study we do not count particles w/o label, because they cannot be scaled
 722           if (!mother)
 723             continue;
 724
 725           // it cannot be just multiplied because we cannot count "half of a particle"
 726           // instead a random generator decides if the particle is counted twice (if value > 1)
 727           // or not (if value < 0)
 728           Int_t bin = fPtSpectrum->FindBin(mother->Pt());
 729           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 730           {
 731             Float_t factor = fPtSpectrum->GetBinContent(bin);
 732             if (factor > 0)
 733             {
 734               Float_t random = gRandom->Uniform();
 735               if (factor > 1 && random < factor - 1)
 736               {
 737                 particleWeight = 2;
 738               }
 739               else if (factor < 1 && random < 1 - factor)
 740                 particleWeight = 0;
 741             }
 742           }
 743         }
 744
 745         //Printf("ESD weight is: %d", particleWeight);
 746
 747         if (TMath::Abs(eta) < 0.5)
 748           nESDTracks05 += particleWeight;
 749
 750         if (TMath::Abs(eta) < 1.0)
 751           nESDTracks10 += particleWeight;
 752
 753         if (TMath::Abs(eta) < 1.3)
 754           nESDTracks14 += particleWeight;
 755
 756         if (fParticleSpecies)
 757         {
 758           Int_t motherLabel = -1;
 759           TParticle* mother = 0;
 760
 761           // find mother
 762           if (label >= 0)
 763             motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
 764           if (motherLabel >= 0)
 765             mother = stack->Particle(motherLabel);
 766
 767           if (!mother)
 768           {
 769             // count tracks that did not have a label
 770             if (TMath::Abs(eta) < etaRange)
 771               nESDTracksSpecies[4]++;
 772           }
 773           else
 774           {
 775             // get particle type (pion, proton, kaon, other) of mother
 776             Int_t idMother = -1;
 777             switch (TMath::Abs(mother->GetPdgCode()))
 778             {
 779               case 211: idMother = 0; break;
 780               case 321: idMother = 1; break;
 781               case 2212: idMother = 2; break;
 782               default: idMother = 3; break;
 783             }
 784
 785             // double counting is ok for particle ratio study
 786             if (TMath::Abs(eta) < etaRange)
 787               nESDTracksSpecies[idMother]++;
 788
 789             // double counting is not ok for efficiency study
 790
 791             // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
 792             if (foundTracks[label])
 793             {
 794               if (TMath::Abs(eta) < etaRange)
 795                 nESDTracksSpecies[6]++;
 796             }
 797             else
 798             {
 799               foundTracks[label] = kTRUE;
 800
 801               // particle (primary) already counted?
 802               if (foundPrimaries[motherLabel])
 803               {
 804                 if (TMath::Abs(eta) < etaRange)
 805                   nESDTracksSpecies[5]++;
 806               }
 807               else
 808                 foundPrimaries[motherLabel] = kTRUE;
 809             }
 810           }
 811         }
 812
 813         if (fParticleCorrection[0])
 814         {
 815           if (label >= 0 && stack->IsPhysicalPrimary(label))
 816           {
 817             TParticle* particle = stack->Particle(label);
 818
 819             // get particle type (pion, proton, kaon, other)
 820             Int_t id = -1;
 821             switch (TMath::Abs(particle->GetPdgCode()))
 822             {
 823               case 211: id = 0; break;
 824               case 321: id = 1; break;
 825               case 2212: id = 2; break;
 826               default: id = 3; break;
 827             }
 828
 829             // todo check if values are not completely off??
 830
 831             // particle (primary) already counted?
 832             if (!foundPrimaries2[label])
 833             {
 834               foundPrimaries2[label] = kTRUE;
 835               fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
 836             }
 837           }
 838         }
 839       }
 840
 841       if (fParticleCorrection[0])
 842       {
 843         // if the particle decays/stops before this radius we do not see it
 844         // 8cm larger than SPD layer 2
 845         // 123cm TPC radius where a track has about 50 clusters (cut limit)
 846         const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
 847
 848         // loop over all primaries that have not been found
 849         for (Int_t i=0; i<nPrim; i++)
 850         {
 851           // already found
 852           if (foundPrimaries2[i])
 853             continue;
 854
 855           TParticle* particle = 0;
 856           TClonesArray* trackrefs = 0;
 857           mcEvent->GetParticleAndTR(i, particle, trackrefs);
 858
 859           // true primary and charged
 860           if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
 861             continue;
 862
 863           //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
 864           if (TMath::Abs(particle->Eta()) > 3)
 865             continue;
 866
 867           // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
 868
 869           // get particle type (pion, proton, kaon, other)
 870           Int_t id = -1;
 871           switch (TMath::Abs(particle->GetPdgCode()))
 872           {
 873             case 211: id = 4; break;
 874             case 321: id = 5; break;
 875             case 2212: id = 6; break;
 876             default: id = 7; break;
 877           }
 878
 879           if (!fParticleCorrection[id])
 880             continue;
 881
 882           // get last track reference
 883           AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
 884
 885           if (!trackref)
 886           {
 887             Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
 888             particle->Print();
 889             continue;
 890           }
 891
 892           // particle in tracking volume long enough...
 893           if (trackref->R() > endRadius)
 894             continue;
 895
 896           if (particle->GetLastDaughter() >= 0)
 897           {
 898             Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
 899             //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
 900             if (uID == kPDecay)
 901             {
 902               // decayed
 903
 904               Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
 905               particle->Print();
 906               Printf("Daugthers:");
 907               for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
 908                 stack->Particle(d)->Print();
 909               Printf(" ");
 910
 911               fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
 912               continue;
 913             }
 914           }
 915
 916           if (trackref->DetectorId() == -1)
 917           {
 918             // stopped
 919             Printf("Particle %d stopped at %f:", i, trackref->R());
 920             particle->Print();
 921             Printf(" ");
 922
 923             fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
 924             continue;
 925           }
 926
 927           Printf("Particle %d simply not tracked", i);
 928           particle->Print();
 929           Printf(" ");
 930         }
 931       }
 932
 933       delete[] foundTracks;
 934       delete[] foundPrimaries;
 935       delete[] foundPrimaries2;
 936
 937 //         if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
 938 //         {
 939 //             TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 940 //             printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
 941 //         }
 942
 943       if (eventTriggered)
 944       {
 945         fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
 946         fMultiplicity->FillNoVertexEvent(vtxMC[2], eventVertex, nMCTracks05, nMCTracks10, nMCTracks14, nESDTracks05, nESDTracks10, nESDTracks14);
 947 //         if (!eventVertex)
 948 //         {
 949 //           if (nESDTracks05 == 0)
 950 //             fTemp1->Fill(nMCTracks05, nESDTracks05);
 951 //         }
 952       }
 953
 954       if (eventTriggered && eventVertex)
 955       {
 956         // fill response matrix using vtxMC (best guess)
 957         fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks14,  nMCTracksAll,  nESDTracks05, nESDTracks10, nESDTracks14);
 958
 959         fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 960
 961         if (fParticleSpecies)
 962           fParticleSpecies->Fill(vtxMC[2], nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
 963       }
 964     }
 965   }
 966
 967   if (etaArr)
 968     delete[] etaArr;
 969   if (labelArr)
 970     delete[] labelArr;
 971
 972 }
 973
 974 void AliMultiplicityTask::Terminate(Option_t *)
 975 {
 976   // The Terminate() function is the last function to be called during
 977   // a query. It always runs on the client, it can be used to present
 978   // the results graphically or save the results to file.
 979
 980   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 981   if (!fOutput) {
 982     Printf("ERROR: fOutput not available");
 983     return;
 984   }
 985
 986   fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
 987   for (Int_t i = 0; i < 8; ++i)
 988     fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
 989   fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
 990
 991   fdNdpT = dynamic_cast<TH1*> (fOutput->FindObject("fdNdpT"));
 992
 993   if (!fMultiplicity)
 994   {
 995     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
 996     return;
 997   }
 998
 999   TString fileName("multiplicity");
1000   if (fSelectProcessType == 1)
1001     fileName += "ND";
1002   if (fSelectProcessType == 2)
1003     fileName += "SD";
1004   if (fSelectProcessType == 3)
1005     fileName += "DD";
1006   fileName += ".root";
1007   TFile* file = TFile::Open(fileName, "RECREATE");
1008
1009   fMultiplicity->SaveHistograms();
1010   for (Int_t i = 0; i < 8; ++i)
1011     if (fParticleCorrection[i])
1012       fParticleCorrection[i]->SaveHistograms();
1013   if (fParticleSpecies)
1014     fParticleSpecies->Write();
1015   if (fdNdpT)
1016     fdNdpT->Write();
1017
1018   fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1019   if (fTemp1)
1020     fTemp1->Write();
1021
1022   fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
1023   if (fEtaPhi)
1024     fEtaPhi->Write();
1025
1026   fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
1027   if (fVertex)
1028     fVertex->Write();
1029
1030   for (Int_t i=0; i<3; i++)
1031   {
1032     fEta[i] = dynamic_cast<TH1*> (fOutput->FindObject(Form("fEta_%d", i)));
1033     if (fEta[i])
1034     {
1035       fEta[i]->Sumw2();
1036       Float_t events = fMultiplicity->GetMultiplicityESD(i)->Integral(1, fMultiplicity->GetMultiplicityESD(i)->GetNbinsX());
1037       if (events > 0)
1038         fEta[i]->Scale(1.0 / events);
1039       fEta[i]->Scale(1.0 / fEta[i]->GetXaxis()->GetBinWidth(1));
1040       fEta[i]->Write();
1041     }
1042   }
1043
1044   TObjString option(fOption);
1045   option.Write();
1046
1047   file->Close();
1048
1049   Printf("Written result to %s", fileName.Data());
1050 }