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 "AliStack.h"
  17 #include "AliESDVertex.h"
  18 #include "AliMultiplicity.h"
  19 #include "AliFMDAnaParameters.h"
  20 //#include "AliFMDParameters.h"
  21
  22 ClassImp(AliFMDAnalysisTaskSharing)
  23
  24 //_____________________________________________________________________
  25 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
  26 : fDebug(0),
  27   fESD(0x0),
  28   foutputESDFMD(),
  29   fEnergy(0),
  30   fNstrips(0),
  31   fSharedThis(kFALSE),
  32   fSharedPrev(kFALSE),
  33   fDiagList(0),
  34   fStandalone(kTRUE),
  35   fEsdVertex(0),
  36   fStatus(kTRUE)
  37 {
  38   // Default constructor
  39   DefineInput (0, AliESDEvent::Class());
  40   DefineOutput(0, AliESDFMD::Class());
  41   DefineOutput(1, AliESDVertex::Class());
  42   DefineOutput(2, AliESDEvent::Class());
  43   DefineOutput(3, TList::Class());
  44 }
  45 //_____________________________________________________________________
  46 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
  47     AliAnalysisTask(name, "AnalysisTaskFMD"),
  48     fDebug(0),
  49     fESD(0x0),
  50     foutputESDFMD(),
  51     fEnergy(0),
  52     fNstrips(0),
  53     fSharedThis(kFALSE),
  54     fSharedPrev(kFALSE),
  55     fDiagList(0),
  56     fStandalone(kTRUE),
  57     fEsdVertex(0),
  58     fStatus(kTRUE)
  59 {
  60   fStandalone = SE;
  61   if(fStandalone) {
  62     DefineInput (0, AliESDEvent::Class());
  63     DefineOutput(0, AliESDFMD::Class());
  64     DefineOutput(1, AliESDVertex::Class());
  65     DefineOutput(2, AliESDEvent::Class());
  66     DefineOutput(3, TList::Class());
  67   }
  68 }
  69 //_____________________________________________________________________
  70 void AliFMDAnalysisTaskSharing::CreateOutputObjects()
  71 {
  72   if(!foutputESDFMD)
  73     foutputESDFMD = new AliESDFMD();
  74
  75   if(!fEsdVertex)
  76     fEsdVertex    = new AliESDVertex();
  77   //Diagnostics
  78   if(!fDiagList)
  79     fDiagList = new TList();
  80
  81   fDiagList->SetName("Sharing diagnostics");
  82   for(Int_t det = 1; det<=3; det++) {
  83     Int_t nRings = (det==1 ? 1 : 2);
  84
  85     for(Int_t iring = 0;iring<nRings; iring++) {
  86       Char_t ringChar = (iring == 0 ? 'I' : 'O');
  87       TH1F* hEdist        = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
  88                                      Form("Edist_before_sharing_FMD%d%c", det, ringChar),
  89                                      1000,0,25);
  90       TH1F* hEdist_after  = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
  91                                      Form("Edist_after_sharing_FMD%d%c", det, ringChar),
  92                                      1000,0,25);
  93
  94
  95       TH1F* hNstripsHit    = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
  96                                      Form("N_strips_hit_FMD%d%c",det,ringChar),
  97                                      25,0,25);
  98       fDiagList->Add(hEdist);
  99       fDiagList->Add(hEdist_after);
 100       fDiagList->Add(hNstripsHit);
 101
 102     }
 103   }
 104 }
 105 //_____________________________________________________________________
 106 void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
 107 {
 108   if(fStandalone)
 109     fESD = (AliESDEvent*)GetInputData(0);
 110 }
 111 //_____________________________________________________________________
 112 void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
 113 {
 114
 115   AliESD* old = fESD->GetAliESDOld();
 116   if (old) {
 117     fESD->CopyFromOldESD();
 118   }
 119
 120   foutputESDFMD->Clear();
 121
 122   Double_t vertex[3];
 123   GetVertex(vertex);
 124   fEsdVertex->SetXYZ(vertex);
 125
 126   if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0) {
 127
 128     fStatus = kFALSE;
 129     return;
 130   }
 131   else
 132     fStatus = kTRUE;
 133   const AliMultiplicity* testmult = fESD->GetMultiplicity();
 134
 135   Int_t nTrackLets = testmult->GetNumberOfTracklets();
 136   if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
 137   else foutputESDFMD->SetUniqueID(kFALSE);
 138
 139   AliESDFMD* fmd = fESD->GetFMDData();
 140   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 141   if (!fmd) return;
 142   Int_t nHits = 0;
 143   for(UShort_t det=1;det<=3;det++) {
 144     Int_t nRings = (det==1 ? 1 : 2);
 145     for (UShort_t ir = 0; ir < nRings; ir++) {
 146       Char_t   ring = (ir == 0 ? 'I' : 'O');
 147       UShort_t nsec = (ir == 0 ? 20  : 40);
 148       UShort_t nstr = (ir == 0 ? 512 : 256);
 149
 150       TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
 151
 152       for(UShort_t sec =0; sec < nsec;  sec++) {
 153         fSharedThis      = kFALSE;
 154         fSharedPrev      = kFALSE;
 155         fEnergy = 0;
 156         fNstrips = 0;
 157         for(UShort_t strip = 0; strip < nstr; strip++) {
 158           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
 159           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 160
 161           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 162
 163           //Double_t eta  = EtaFromStrip(det,ring,sec,strip,vertex[2]);//fmd->Eta(det,ring,sec,strip);
 164           //Double_t eta = fmd->Eta(det,ring,sec,strip);
 165           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 166           //std::cout<<EtaFromStrip(det,ring,sec,strip,vertex[2]) <<"    "<<fmd->Eta(det,ring,sec,strip)<<std::endl;
 167
 168           hEdist->Fill(mult);
 169           if(fmd->IsAngleCorrected())
 170             mult = mult/TMath::Cos(Eta2Theta(eta));
 171           Float_t Eprev = 0;
 172           Float_t Enext = 0;
 173           if(strip != 0)
 174             if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
 175               Eprev = fmd->Multiplicity(det,ring,sec,strip-1);
 176               if(fmd->IsAngleCorrected())
 177                 Eprev = Eprev/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
 178             }
 179           if(strip != nstr - 1)
 180             if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
 181               Enext = fmd->Multiplicity(det,ring,sec,strip+1);
 182               if(fmd->IsAngleCorrected())
 183                 Enext = Enext/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
 184             }
 185
 186           Float_t merged_energy = GetMultiplicityOfStrip(mult,eta,Eprev,Enext,det,ring,sec,strip);
 187
 188           if(merged_energy > 0 )
 189             nHits++;
 190           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,merged_energy);
 191           foutputESDFMD->SetEta(det,ring,sec,strip,eta);
 192
 193         }
 194       }
 195     }
 196   }
 197
 198   if(fStandalone) {
 199     PostData(0, foutputESDFMD);
 200     PostData(1, fEsdVertex);
 201     PostData(2, fESD);
 202     PostData(3, fDiagList);
 203   }
 204 }
 205 //_____________________________________________________________________
 206 Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
 207                                                           Float_t eta,
 208                                                           Float_t Eprev,
 209                                                           Float_t Enext,
 210                                                           UShort_t   det,
 211                                                           Char_t  ring,
 212                                                           UShort_t /*sec*/,
 213                                                           UShort_t /*strip*/) {
 214   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 215
 216   Float_t merged_energy = 0;
 217   //Float_t nParticles = 0;
 218   Float_t cutLow  = 0.15;
 219   if(ring == 'I')
 220     cutLow = 0.1;
 221
 222   //cutLow = 0;
 223   //AliFMDParameters* recopars = AliFMDParameters::Instance();
 224   //cutLow = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
 225
 226
 227
 228   Float_t cutHigh = pars->GetMPV(det,ring,eta) - 3*pars->GetSigma(det,ring,eta);
 229   // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
 230   Float_t Etotal  = mult;
 231
 232   //if(mult > 5)
 233   //  std::cout<<mult<<"    "<<det<<"    "<<ring<<"   "<<sec<<"    "<<strip<<std::endl;
 234   //Float_t slow_particle_cut = 2*pars->GetMPV(det,ring,eta);
 235
 236   //if(recopars->IsDead(det,ring,sec,strip))
 237   //  std::cout<<"dead channel"<<std::endl;
 238   //if(foutputESDFMD->GetUniqueID() == kTRUE) {
 239   // Float_t mpv = pars->GetMPV(det,ring,eta);
 240
 241   /*
 242   if(foutputESDFMD->GetUniqueID() == kFALSE) {
 243
 244     if(mult > 15)
 245       return 0;
 246
 247     if(mult > cutLow ) {
 248       fEnergy = fEnergy + mult;
 249       fNstrips++;
 250     }
 251
 252     if( (Enext < cutLow && fEnergy > 0 )   || fNstrips >2 ){
 253
 254       //if((fEnergy*TMath::Cos(Eta2Theta(eta))) > cutPart || fNstrips > 1) {
 255       nParticles = 1;
 256       merged_energy = fEnergy*TMath::Cos(Eta2Theta(eta));
 257       TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
 258       hEdist->Fill(fEnergy);
 259       TH1F* hNstrips = (TH1F*)fDiagList.FindObject(Form("N_strips_hit_FMD%d%c",det,ring));
 260       hNstrips->Fill(fNstrips);
 261       // std::cout<<Form("Merged signals %f %f %f into %f , %f in det %d, ring %c, sec %d, strip %d",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,det,ring,sec,strip )<<std::endl;
 262
 263       // }
 264     // else
 265     //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, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
 266
 267       fEnergy  = 0;
 268       fNstrips = 0;
 269
 270     return merged_energy;
 271     }
 272
 273
 274     return 0;
 275
 276   }
 277   else {*/
 278     //std::cout<<det<<ring<<"   "<<sec<<"    "<<strip<<"   "<<cutLow<<std::endl;
 279   if(fSharedThis) {
 280     fSharedThis      = kFALSE;
 281     fSharedPrev      = kTRUE;
 282     return 0.;
 283   }
 284
 285   /*  if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
 286     fSharedThis      = kFALSE;
 287     fSharedPrev      = kFALSE;
 288     return 0;
 289     }*/
 290   if(mult<Enext && Enext>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
 291     {
 292       fSharedThis      = kFALSE;
 293       fSharedPrev      = kFALSE;
 294       return 0;
 295     }
 296   if(mult > 15)
 297     {
 298       //   std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
 299       fSharedThis      = kFALSE;
 300       fSharedPrev      = kFALSE;
 301       return 0;
 302     }
 303
 304   if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
 305     Etotal += Eprev;
 306   }
 307
 308   if(Enext > cutLow && Enext < cutHigh ) {
 309     Etotal += Enext;
 310     fSharedThis      = kTRUE;
 311   }
 312   TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
 313   hEdist->Fill(Etotal);
 314
 315   Etotal = Etotal*TMath::Cos(Eta2Theta(eta));
 316   if(Etotal > 0) {
 317
 318     merged_energy = Etotal;
 319     fSharedPrev      = kTRUE;
 320     // if(det == 1 && ring =='I')
 321     // 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;
 322   }
 323     else{// if(Etotal > 0) {
 324       //if(det == 3 && ring =='I')
 325       //        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;
 326     fSharedThis      = kFALSE;
 327     fSharedPrev      = kFALSE;
 328   }
 329   // merged_energy = mult;
 330
 331   return merged_energy;
 332   //}
 333 }
 334 //_____________________________________________________________________
 335 void AliFMDAnalysisTaskSharing::GetVertex(Double_t* vertexXYZ)
 336 {
 337   const AliESDVertex* vertex = 0;
 338   vertex = fESD->GetPrimaryVertex();
 339   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 340     vertex = fESD->GetPrimaryVertexSPD();
 341   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 342     vertex = fESD->GetPrimaryVertexTPC();
 343   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 344     vertex = fESD->GetVertex();
 345   if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
 346     vertex->GetXYZ(vertexXYZ);
 347     //std::cout<<vertex->GetName()<<"   "<< vertex->GetTitle() <<"   "<< vertex->GetZv()<<std::endl;
 348     return;
 349   }
 350   else if (fESD->GetESDTZERO()) {
 351     vertexXYZ[0] = 0;
 352     vertexXYZ[1] = 0;
 353     vertexXYZ[2] = fESD->GetT0zVertex();
 354
 355     return;
 356   }
 357
 358   return;
 359
 360 }
 361 //_____________________________________________________________________
 362 Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) {
 363
 364   Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
 365
 366   if(eta < 0)
 367     theta = theta-TMath::Pi();
 368
 369   //  std::cout<<"From eta2Theta: "<<theta<<"   "<<eta<<std::endl;
 370   return theta;
 371
 372
 373
 374 }
 375 //_____________________________________________________________________
 376 /*Double_t AliFMDAnalysisTaskSharing::EtaFromStrip(UShort_t det,
 377                                                 Char_t ring,
 378                                                 UShort_t sector,
 379                                                 UShort_t strip,
 380                                                 Double_t zvtx)
 381 {
 382
 383   AliFMDGeometry* geo = AliFMDGeometry::Instance();
 384
 385   Double_t x,y,z;
 386   geo->Detector2XYZ(det,ring,sector,strip,x,y,z);
 387
 388   Double_t r = TMath::Sqrt(x*x+y*y);
 389
 390   Double_t z_real      = z-zvtx;
 391   Double_t theta       = TMath::ATan2(r,z_real);
 392   // std::cout<<"From EtaFromStrip "<<theta<<std::endl;
 393   Double_t eta         =  -1*TMath::Log(TMath::Tan(0.5*theta));
 394
 395   // std::cout<<det<<"   "<<ring<<"   "<<sector<<"   "<<strip<<"   "<<r<<"    "<<z_real<<"   "<<theta<<"    "<<eta<<std::endl;
 396
 397   return eta;
 398   }*/
 399 //_____________________________________________________________________
 400 //
 401 // EOF
 402 //