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