PWG2/FORWARD/analysis/AliFMDAnalysisTaskSharing.cxx

   1
   2 #include <TROOT.h>
   3 #include <TSystem.h>
   4 #include <TInterpreter.h>
   5 #include <TChain.h>
   6 #include <TFile.h>
   7 #include <TList.h>
   8 #include <iostream>
   9 #include <TMath.h>
  10 //#include "AliFMDDebug.h"
  11 #include "AliFMDAnalysisTaskSharing.h"
  12 #include "AliAnalysisManager.h"
  13 #include "AliESDFMD.h"
  14 //#include "AliFMDGeometry.h"
  15 #include "AliMCEventHandler.h"
  16 #include "AliMCEvent.h"
  17 #include "AliStack.h"
  18 #include "AliESDVertex.h"
  19 #include "AliMultiplicity.h"
  20 #include "AliFMDAnaParameters.h"
  21 #include "TH1F.h"
  22 #include "TObjString.h"
  23 //#include "/home/canute/ALICE/AliRoot/PWG0/AliPWG0Helper.h"
  24 //#include "AliFMDParameters.h"
  25 #include "AliGenEventHeader.h"
  26 #include "AliHeader.h"
  27 #include "AliStack.h"
  28 #include "AliMCParticle.h"
  29 #include "AliFMDStripIndex.h"
  30
  31 // This is the task to do the FMD sharing or hit merging.
  32 // It reads the input ESDFMD data and posts an ESDFMD object to
  33 // the tasks that must be performed after this task ie.
  34 // Density, BackgroundCorrection and Dndeta.
  35 // Author: Hans Hjersing Dalsgaard, hans.dalsgaard@cern.ch
  36
  37
  38 ClassImp(AliFMDAnalysisTaskSharing)
  39
  40 //_____________________________________________________________________
  41 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
  42 : fDebug(0),
  43   fESD(0x0),
  44   foutputESDFMD(),
  45   fSharedThis(kFALSE),
  46   fSharedPrev(kFALSE),
  47   fDiagList(0),
  48   fStandalone(kTRUE),
  49   fEsdVertex(0),
  50   fStatus(kTRUE),
  51   fLastTrackByStrip(0)
  52 {
  53   // Default constructor
  54   DefineInput (0, AliESDEvent::Class());
  55   DefineOutput(0, AliESDFMD::Class());
  56   DefineOutput(1, AliESDVertex::Class());
  57   DefineOutput(2, AliESDEvent::Class());
  58   DefineOutput(3, TList::Class());
  59 }
  60 //_____________________________________________________________________
  61 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
  62     AliAnalysisTask(name, "AnalysisTaskFMD"),
  63     fDebug(0),
  64     fESD(0x0),
  65     foutputESDFMD(),
  66     fSharedThis(kFALSE),
  67     fSharedPrev(kFALSE),
  68     fDiagList(0),
  69     fStandalone(kTRUE),
  70     fEsdVertex(0),
  71     fStatus(kTRUE),
  72     fLastTrackByStrip(0)
  73 {
  74   // named constructor
  75   fStandalone = SE;
  76   if(fStandalone) {
  77     DefineInput (0, AliESDEvent::Class());
  78     DefineOutput(0, AliESDFMD::Class());
  79     DefineOutput(1, AliESDVertex::Class());
  80     DefineOutput(2, AliESDEvent::Class());
  81     DefineOutput(3, TList::Class());
  82   }
  83 }
  84 //_____________________________________________________________________
  85 void AliFMDAnalysisTaskSharing::CreateOutputObjects()
  86 {
  87   // Create the output objects
  88   if(!foutputESDFMD)
  89     foutputESDFMD = new AliESDFMD();
  90
  91   if(!fEsdVertex)
  92     fEsdVertex    = new AliESDVertex();
  93   //Diagnostics
  94   if(!fDiagList)
  95     fDiagList = new TList();
  96
  97   fDiagList->SetName("Sharing diagnostics");
  98
  99   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 100   TH2F* hBg   = pars->GetBackgroundCorrection(1, 'I', 0);
 101   TH1F* hPrimary = new TH1F("hMultvsEtaNoCuts","hMultvsEtaNoCuts",
 102                             hBg->GetNbinsX(),
 103                             hBg->GetXaxis()->GetXmin(),
 104                             hBg->GetXaxis()->GetXmax());
 105   hPrimary->Sumw2();
 106   fDiagList->Add(hPrimary);
 107   TH1F* hZvtx = new TH1F("hZvtx","z vertex distribution",pars->GetNvtxBins(),-1*pars->GetVtxCutZ(),pars->GetVtxCutZ());
 108
 109   fDiagList->Add(hZvtx);
 110
 111   TH1F* hPrimVertexBin = 0;
 112   TH1F* hHits = 0;
 113   for(Int_t i = 0; i< pars->GetNvtxBins(); i++) {
 114
 115     hPrimVertexBin = new TH1F(Form("primmult_NoCuts_vtxbin%d",i),
 116                               Form("primmult_NoCuts_vtxbin%d",i),
 117                               hBg->GetNbinsX(),
 118                               hBg->GetXaxis()->GetXmin(),
 119                               hBg->GetXaxis()->GetXmax());
 120     hPrimVertexBin->Sumw2();
 121     fDiagList->Add(hPrimVertexBin);
 122
 123   }
 124
 125   for(Int_t det = 1; det<=3; det++) {
 126     Int_t nRings = (det==1 ? 1 : 2);
 127
 128     for(Int_t iring = 0;iring<nRings; iring++) {
 129       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 130       TH1F* hEdist        = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
 131                                      Form("Edist_before_sharing_FMD%d%c", det, ringChar),
 132                                      1000,0,25);
 133       TH1F* hEdistAfter  = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
 134                                      Form("Edist_after_sharing_FMD%d%c", det, ringChar),
 135                                      1000,0,25);
 136
 137
 138       //TH1F* hNstripsHit    = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
 139       //                                     Form("N_strips_hit_FMD%d%c",det,ringChar),
 140       //                                     25,0,25);
 141       fDiagList->Add(hEdist);
 142       fDiagList->Add(hEdistAfter);
 143       //fDiagList->Add(hNstripsHit);
 144
 145       for(Int_t i = 0; i< pars->GetNvtxBins(); i++) {
 146         hHits  = new TH1F(Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
 147                           hBg->GetNbinsX(),
 148                           hBg->GetXaxis()->GetXmin(),
 149                           hBg->GetXaxis()->GetXmax());
 150         hHits->Sumw2();
 151         fDiagList->Add(hHits);
 152
 153       }
 154
 155     }
 156   }
 157   TH1F*  nMCevents = new TH1F("nMCEventsNoCuts","nMCEventsNoCuts",pars->GetNvtxBins(),0,pars->GetNvtxBins());
 158
 159   fDiagList->Add(nMCevents);
 160
 161 }
 162 //_____________________________________________________________________
 163 void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
 164 {
 165   // connect the input data
 166   if(fStandalone)
 167     fESD = (AliESDEvent*)GetInputData(0);
 168 }
 169 //_____________________________________________________________________
 170 void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
 171 {
 172   //perform analysis on one event
 173   AliESD* old = fESD->GetAliESDOld();
 174   if (old) {
 175     fESD->CopyFromOldESD();
 176   }
 177
 178   foutputESDFMD->Clear();
 179
 180   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 181   Double_t vertex[3];
 182   Bool_t vtxStatus = pars->GetVertex(fESD,vertex);
 183   fEsdVertex->SetXYZ(vertex);
 184
 185   // Process primaries here to get true MC distribution
 186   if(pars->GetProcessPrimary())
 187     ProcessPrimary();
 188
 189   Bool_t isTriggered = pars->IsEventTriggered(fESD);
 190
 191   if(!isTriggered) {
 192     fStatus = kFALSE;
 193     return;
 194    }
 195    else
 196      fStatus = kTRUE;
 197
 198   if(!vtxStatus) {
 199     fStatus = kFALSE;
 200     return;
 201   }
 202   else
 203     fStatus = kTRUE;
 204
 205   TH1F* hZvtx = (TH1F*)fDiagList->FindObject("hZvtx");
 206   hZvtx->Fill(vertex[2]);
 207
 208
 209   const AliMultiplicity* testmult = fESD->GetMultiplicity();
 210
 211   Int_t nTrackLets = testmult->GetNumberOfTracklets();
 212   if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
 213   else foutputESDFMD->SetUniqueID(kFALSE);
 214
 215   AliESDFMD* fmd = fESD->GetFMDData();
 216
 217   if (!fmd) return;
 218   Int_t nHits = 0;
 219   for(UShort_t det=1;det<=3;det++) {
 220     Int_t nRings = (det==1 ? 1 : 2);
 221     for (UShort_t ir = 0; ir < nRings; ir++) {
 222       Char_t   ring = (ir == 0 ? 'I' : 'O');
 223       UShort_t nsec = (ir == 0 ? 20  : 40);
 224       UShort_t nstr = (ir == 0 ? 512 : 256);
 225
 226       TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
 227
 228       for(UShort_t sec =0; sec < nsec;  sec++) {
 229         fSharedThis      = kFALSE;
 230         fSharedPrev      = kFALSE;
 231
 232         for(UShort_t strip = 0; strip < nstr; strip++) {
 233           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
 234           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 235
 236           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 237
 238           //Double_t eta  = EtaFromStrip(det,ring,sec,strip,vertex[2]);//fmd->Eta(det,ring,sec,strip);
 239           //Double_t eta = fmd->Eta(det,ring,sec,strip);
 240           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 241           //std::cout<<EtaFromStrip(det,ring,sec,strip,vertex[2]) <<"    "<<fmd->Eta(det,ring,sec,strip)<<std::endl;
 242
 243           hEdist->Fill(mult);
 244           if(fmd->IsAngleCorrected())
 245             mult = mult/TMath::Cos(Eta2Theta(eta));
 246           Float_t prevE = 0;
 247           Float_t nextE = 0;
 248           if(strip != 0)
 249             if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
 250               prevE = fmd->Multiplicity(det,ring,sec,strip-1);
 251               if(fmd->IsAngleCorrected())
 252                 prevE = prevE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
 253             }
 254           if(strip != nstr - 1)
 255             if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
 256               nextE = fmd->Multiplicity(det,ring,sec,strip+1);
 257               if(fmd->IsAngleCorrected())
 258                 nextE = nextE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
 259             }
 260
 261           Float_t mergedEnergy = GetMultiplicityOfStrip(mult,eta,prevE,nextE,det,ring,sec,strip);
 262
 263           if(mergedEnergy > 0 )
 264             nHits++;
 265           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,mergedEnergy);
 266           foutputESDFMD->SetEta(det,ring,sec,strip,eta);
 267
 268         }
 269       }
 270     }
 271   }
 272
 273
 274   if(fStandalone) {
 275     PostData(0, foutputESDFMD);
 276     PostData(1, fEsdVertex);
 277     PostData(2, fESD);
 278     PostData(3, fDiagList);
 279   }
 280 }
 281 //_____________________________________________________________________
 282 Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
 283                                                           Float_t eta,
 284                                                           Float_t prevE,
 285                                                           Float_t nextE,
 286                                                           UShort_t   det,
 287                                                           Char_t  ring,
 288                                                           UShort_t /*sec*/,
 289                                                           UShort_t /*strip*/) {
 290   //analyse and perform sharing on one strip
 291   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 292
 293   Float_t mergedEnergy = 0;
 294   //Float_t nParticles = 0;
 295   Float_t cutLow  = 0.25;//0.15;
 296   // if(ring == 'I')
 297   //  cutLow = 0.1;
 298
 299   //cutLow = 0;
 300   //AliFMDParameters* recopars = AliFMDParameters::Instance();
 301   //cutLow = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
 302
 303
 304
 305   Float_t cutHigh = pars->GetMPV(det,ring,eta) - 3*pars->GetSigma(det,ring,eta);
 306
 307   // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
 308   Float_t totalE  = mult;
 309
 310
 311     //std::cout<<det<<ring<<"   "<<sec<<"    "<<strip<<"   "<<cutLow<<std::endl;
 312   if(fSharedThis) {
 313     fSharedThis      = kFALSE;
 314     fSharedPrev      = kTRUE;
 315     return 0.;
 316   }
 317
 318   /*  if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
 319     fSharedThis      = kFALSE;
 320     fSharedPrev      = kFALSE;
 321     return 0;
 322     }*/
 323   if(mult<nextE && nextE>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
 324     {
 325       fSharedThis      = kFALSE;
 326       fSharedPrev      = kFALSE;
 327       return 0;
 328     }
 329   if(mult > 15)
 330     {
 331       //   std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
 332       fSharedThis      = kFALSE;
 333       fSharedPrev      = kFALSE;
 334       return 0;
 335     }
 336
 337   if(prevE > cutLow && prevE < cutHigh && !fSharedPrev ) {
 338     totalE += prevE;
 339   }
 340
 341   if(nextE > cutLow && nextE < cutHigh ) {
 342     totalE += nextE;
 343     fSharedThis      = kTRUE;
 344   }
 345   TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
 346   hEdist->Fill(totalE);
 347
 348   totalE = totalE*TMath::Cos(Eta2Theta(eta));
 349   if(totalE > 0) {
 350
 351     mergedEnergy = totalE;
 352     fSharedPrev      = kTRUE;
 353     // if(det == 1 && ring =='I')
 354     // std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",prevE, mult, nextE, totalE/TMath::Cos(Eta2Theta(eta)),totalE,strip,sec,ring,det )<<std::endl;
 355   }
 356     else{// if(totalE > 0) {
 357       //if(det == 3 && ring =='I')
 358       //        std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",prevE, mult, nextE, totalE/TMath::Cos(Eta2Theta(eta)),totalE,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
 359     fSharedThis      = kFALSE;
 360     fSharedPrev      = kFALSE;
 361   }
 362   // mergedEnergy = mult;
 363
 364   return mergedEnergy;
 365   //}
 366 }
 367
 368 //_____________________________________________________________________
 369 Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) const{
 370   //convert the eta of a strip to a theta
 371   Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
 372
 373   if(eta < 0)
 374     theta = theta-TMath::Pi();
 375
 376   //  std::cout<<"From eta2Theta: "<<theta<<"   "<<eta<<std::endl;
 377   return theta;
 378
 379
 380
 381 }
 382
 383
 384
 385 //_____________________________________________________________________
 386 void AliFMDAnalysisTaskSharing::ProcessPrimary() {
 387   //Get the undspoiled MC dN/deta before event cuts
 388   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 389   AliMCEvent* mcEvent = eventHandler->MCEvent();
 390   if(!mcEvent)
 391     return;
 392   fLastTrackByStrip.Reset(-1);
 393   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 394
 395   AliMCParticle* particle = 0;
 396
 397   AliStack* stack = mcEvent->Stack();
 398
 399   TH1F* hPrimary = (TH1F*)fDiagList->FindObject("hMultvsEtaNoCuts");
 400   AliHeader* header            = mcEvent->Header();
 401   AliGenEventHeader* genHeader = header->GenEventHeader();
 402
 403
 404
 405   TArrayF vertex;
 406   genHeader->PrimaryVertex(vertex);
 407
 408   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
 409     return;
 410
 411   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 412   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
 413   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 414
 415   Bool_t firstTrack = kTRUE;
 416
 417   Int_t nTracks = stack->GetNprimary();
 418   if(pars->GetProcessHits())
 419     nTracks = stack->GetNtrack();
 420   TH1F* nMCevents = (TH1F*)fDiagList->FindObject("nMCEventsNoCuts");
 421   for(Int_t i = 0 ;i<nTracks;i++) {
 422     particle = (AliMCParticle*) mcEvent->GetTrack(i);
 423     if(!particle)
 424       continue;
 425
 426     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 427       hPrimary->Fill(particle->Eta());
 428
 429
 430       TH1F* hPrimVtxBin = (TH1F*)fDiagList->FindObject(Form("primmult_NoCuts_vtxbin%d",vertexBin));
 431       hPrimVtxBin->Fill(particle->Eta());
 432
 433       if(firstTrack) {
 434         nMCevents->Fill(vertexBin);
 435         firstTrack = kFALSE;
 436       }
 437
 438     }
 439      if(pars->GetProcessHits()) {
 440
 441       for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 442
 443         AliTrackReference* ref = particle->GetTrackReference(j);
 444         UShort_t det,sec,strip;
 445         Char_t   ring;
 446         if(ref->DetectorId() != AliTrackReference::kFMD)
 447           continue;
 448         AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 449         Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
 450         if(particle->Charge() != 0 && i != thisStripTrack ) {
 451           //Double_t x,y,z;
 452
 453           Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
 454           TH1F* hHits = (TH1F*)fDiagList->FindObject(Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vertexBin));
 455
 456
 457           hHits->Fill(eta);
 458
 459           Float_t nstrips = (ring =='O' ? 256 : 512);
 460
 461           fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
 462
 463           if(strip >0)
 464             fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
 465           if(strip < (nstrips - 1))
 466             fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
 467
 468         }
 469
 470
 471       }
 472
 473
 474     }
 475
 476   }
 477
 478 }
 479
 480 //_____________________________________________________________________
 481 //
 482 // EOF
 483 //