]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MFT/AliMFT.cxx
Fix for ESD analysis
[u/mrichter/AliRoot.git] / MFT / AliMFT.cxx
index 1a23486c8eef49d845ed907f5c21743a4e8a1c3b..f20a504454517d82f65cc2f3c85c192b05cd1493 100644 (file)
@@ -53,17 +53,17 @@ AliMFT::AliMFT():
   AliDetector(),
   fVersion(1),
   fNPlanes(0),
-  fNSlices(0),
+  fNSlices(1),
   fSDigitsPerPlane(0),
   fDigitsPerPlane(0),
   fRecPointsPerPlane(0),
   fSideDigits(0),
   fSegmentation(0),
   fNameGeomFile(0),
-  fChargeDispersion(0),
+  fChargeDispersion(25.e-4),
   fSingleStepForChargeDispersion(0),
-  fNStepForChargeDispersion(0),
-  fDensitySiOverSupport(10)
+  fNStepForChargeDispersion(4),
+  fDensitySupportOverSi(0.15)
 {
 
   // default constructor
@@ -76,17 +76,17 @@ AliMFT::AliMFT(const Char_t *name, const Char_t *title):
   AliDetector(name, title),
   fVersion(1),
   fNPlanes(0),
-  fNSlices(0),
+  fNSlices(1),
   fSDigitsPerPlane(0),
   fDigitsPerPlane(0),
   fRecPointsPerPlane(0),
   fSideDigits(0),
   fSegmentation(0),
   fNameGeomFile(0),
-  fChargeDispersion(0),
+  fChargeDispersion(25.e-4),
   fSingleStepForChargeDispersion(0),
-  fNStepForChargeDispersion(0),
-  fDensitySiOverSupport(10)
+  fNStepForChargeDispersion(4),
+  fDensitySupportOverSi(0.15)
 {
 
   fNameGeomFile = "AliMFTGeometry.root";
@@ -103,17 +103,17 @@ AliMFT::AliMFT(const Char_t *name, const Char_t *title, Char_t *nameGeomFile):
   AliDetector(name, title),
   fVersion(1),
   fNPlanes(0),
-  fNSlices(0),
+  fNSlices(1),
   fSDigitsPerPlane(0),
   fDigitsPerPlane(0),
   fRecPointsPerPlane(0),
   fSideDigits(0),
   fSegmentation(0),
   fNameGeomFile(0),
-  fChargeDispersion(0),
+  fChargeDispersion(25.e-4),
   fSingleStepForChargeDispersion(0),
-  fNStepForChargeDispersion(0),
-  fDensitySiOverSupport(10)
+  fNStepForChargeDispersion(4),
+  fDensitySupportOverSi(0.15)
 {
 
   fNameGeomFile = nameGeomFile;
@@ -175,7 +175,7 @@ void AliMFT::CreateMaterials() {
   AliMaterial(++matId, "Readout", aSi,   zSi,    dSi,    radSi,  absSi  );  
   AliMedium(kReadout,  "Readout", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
 
-  AliMaterial(++matId, "Support", aSi,   zSi,    dSi/fDensitySiOverSupport, fDensitySiOverSupport*radSi, fDensitySiOverSupport*absSi);  
+  AliMaterial(++matId, "Support", aSi,   zSi,    dSi*fDensitySupportOverSi, radSi/fDensitySupportOverSi, absSi/fDensitySupportOverSi);  
   AliMedium(kSupport,  "Support", matId, unsens, isxfld,  sxmgmx,    tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
     
   AliInfo("End MFT materials");
@@ -201,10 +201,29 @@ void AliMFT::CreateGeometry() {
 
 //====================================================================================================================================================
 
+void AliMFT::AddAlignableVolumes() {
+
+  // Create entries for alignable volumes associating the symbolic volume
+  // name with the corresponding volume path. Needs to be syncronized with
+  // eventual changes in the geometry.
+
+  TString sysName = "MFT";
+  TString volPath = "/ALIC_1/MFT_0";
+  
+  if (!gGeoManager->SetAlignableEntry(sysName.Data(),volPath.Data())) {
+    AliFatal(Form("Alignable entry %s not created. Volume path %s not valid", sysName.Data(), volPath.Data()));
+  }  
+
+}
+
+//====================================================================================================================================================
+
 void AliMFT::StepManager() {
 
   // Full Step Manager
 
+  AliDebug(2, Form("Entering StepManager: gMC->CurrentVolName() = %s", gMC->CurrentVolName()));
+
   if (!fSegmentation) AliFatal("No segmentation available");    // DO WE HAVE A SEGMENTATION???
 
   if (!(this->IsActive())) return;
@@ -247,8 +266,8 @@ void AliMFT::StepManager() {
   gMC->TrackPosition(position);
   gMC->TrackMomentum(momentum);
 
-  AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (z=%f) belongs to track %02d\n", 
-                  gMC->CurrentVolName(), fNhits, position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber())); 
+  AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (x=%f, y=%f, z=%f) belongs to track %02d\n", 
+                  gMC->CurrentVolName(), fNhits, position.X(), position.Y(), position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber())); 
 
   hit.SetPosition(position);
   hit.SetTOF(gMC->TrackTime());
@@ -309,6 +328,8 @@ TGeoVolumeAssembly* AliMFT::CreateVol() {
                                                    plane->GetRMaxSupport(),
                                                    0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() - 
                                                         plane->GetSupportElement(0)->GetAxis(2)->GetXmin()) );
+    AliDebug(2, Form("Created vol %s", supportElem->GetName()));
+    supportElem->SetLineColor(kCyan-9);
     vol -> AddNode(supportElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
 
     AliDebug(1, "support elements created!");
@@ -328,6 +349,8 @@ TGeoVolumeAssembly* AliMFT::CreateVol() {
       for (Int_t iSlice=0; iSlice<fNSlices; iSlice++) {
        origin[2] = -0.5*(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice+1) + plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice) );
        TGeoVolume *activeElem = gGeoManager->MakeBox(Form("MFT_plane%02d_active%03d_slice%02d", iPlane, iActive, iSlice), silicon, dx, dy, dz);
+       AliDebug(2, Form("Created vol %s", activeElem->GetName()));
+       activeElem->SetLineColor(kGreen);
        vol -> AddNode(activeElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
       }
 
@@ -348,6 +371,8 @@ TGeoVolumeAssembly* AliMFT::CreateVol() {
       origin[2] = -0.5*(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
 
       TGeoVolume *readoutElem = gGeoManager->MakeBox(Form("MFT_plane%02d_readout%03d", iPlane, iReadout), readout, dx, dy, dz);
+      AliDebug(2, Form("Created vol %s", readoutElem->GetName()));
+      readoutElem->SetLineColor(kRed);
       vol -> AddNode(readoutElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
       
     }
@@ -408,7 +433,7 @@ void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t
 
   //  Add sdigits of these hits to the list
   
-  AliDebug(1, "Start Hits2SDigitsLocal");
+  AliDebug(1, "Entering Hits2SDigitsLocal");
   
   if (!fSegmentation) CreateGeometry();
   
@@ -424,82 +449,155 @@ void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t
 
     AliMFTHit *hit = (AliMFTHit*) hits->At(iHit);
 
-    AliMFTDigit sDigit;
-    sDigit.SetEloss(hit->GetEloss());
-    sDigit.SetDetElemID(hit->GetDetElemID());
-    sDigit.SetPlane(hit->GetPlane());
-    sDigit.AddMCLabel(hit->GetTrack()); 
+    // Creating "main digit"
+
+    AliMFTDigit *mainSDigit = new AliMFTDigit();
+    mainSDigit->SetEloss(hit->GetEloss());
+    mainSDigit->SetDetElemID(hit->GetDetElemID());
+    mainSDigit->SetPlane(hit->GetPlane());
+    mainSDigit->AddMCLabel(hit->GetTrack()); 
 
     Int_t xPixel = -1;
     Int_t yPixel = -1;
-    if (fSegmentation->Hit2PixelID(hit->X(), hit->Y(), sDigit.GetDetElemID(), xPixel, yPixel)) {
-      sDigit.SetPixID(xPixel, yPixel, 0);
-      sDigit.SetPixWidth(fSegmentation->GetPixelSizeX(sDigit.GetDetElemID()), 
-                        fSegmentation->GetPixelSizeY(sDigit.GetDetElemID()),
-                        fSegmentation->GetPixelSizeZ(sDigit.GetDetElemID()));  
-      sDigit.SetPixCenter(fSegmentation->GetPixelCenterX(sDigit.GetDetElemID(), xPixel), 
-                         fSegmentation->GetPixelCenterY(sDigit.GetDetElemID(), yPixel),
-                         fSegmentation->GetPixelCenterZ(sDigit.GetDetElemID(), 0));  
-      new ((*pSDigList[sDigit.GetPlane()])[pSDigList[sDigit.GetPlane()]->GetEntries()]) AliMFTDigit(sDigit);
+    if (fSegmentation->Hit2PixelID(hit->X(), hit->Y(), mainSDigit->GetDetElemID(), xPixel, yPixel)) {
+      mainSDigit->SetPixID(xPixel, yPixel, 0);
+      mainSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(mainSDigit->GetDetElemID()), 
+                         fSegmentation->GetPixelSizeY(mainSDigit->GetDetElemID()),
+                         fSegmentation->GetPixelSizeZ(mainSDigit->GetDetElemID()));  
+      mainSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(mainSDigit->GetDetElemID(), xPixel), 
+                          fSegmentation->GetPixelCenterY(mainSDigit->GetDetElemID(), yPixel),
+                          fSegmentation->GetPixelCenterZ(mainSDigit->GetDetElemID(), 0));
+      new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*mainSDigit);
       AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
-                      sDigit.GetPixelCenterX(), sDigit.GetPixelCenterY(), sDigit.GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
-//       AliDebug(1, Form("Created new sdigit from hit: residual is (%f, %f, %f)",
-//                    sDigit.GetPixelCenterX()-hit->X(), sDigit.GetPixelCenterY()-hit->Y(), sDigit.GetPixelCenterZ()-hit->Z()));
+                              mainSDigit->GetPixelCenterX(), mainSDigit->GetPixelCenterY(), mainSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
     }
 
-    // creating "side hits" to simulate the effect of charge dispersion
+    // creating "side digits" to simulate the effect of charge dispersion
 
-    Int_t xPixelNew = -1;
-    Int_t yPixelNew = -1;
-    Double_t x0 = hit->X();
-    Double_t y0 = hit->Y();
     Double_t pi4 = TMath::Pi()/4.;
     for (Int_t iStep=0; iStep<fNStepForChargeDispersion; iStep++) {
       Double_t shift = (iStep+1) * fSingleStepForChargeDispersion;
       for (Int_t iAngle=0; iAngle<8; iAngle++) {
        Double_t shiftX = shift*TMath::Cos(iAngle*pi4);
        Double_t shiftY = shift*TMath::Sin(iAngle*pi4);
-       if (fSegmentation->Hit2PixelID(x0+shiftX, y0+shiftY, hit->GetDetElemID(), xPixelNew, yPixelNew)) {
+       if (fSegmentation->Hit2PixelID(hit->X()+shiftX, hit->Y()+shiftY, hit->GetDetElemID(), xPixel, yPixel)) {
          Bool_t digitExists = kFALSE;
-         if (xPixelNew==xPixel && yPixelNew==yPixel) digitExists = kTRUE;
-         if (!digitExists) {
-           for (Int_t iSideDigit=0; iSideDigit<fSideDigits->GetEntries(); iSideDigit++) {
-             if (xPixelNew==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelX() && 
-                 yPixelNew==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelY()) {
-               digitExists = kTRUE;
-               break;
-             }
+         for (Int_t iSideDigit=0; iSideDigit<fSideDigits->GetEntries(); iSideDigit++) {
+           if (xPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelX() && 
+               yPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelY()) {
+             digitExists = kTRUE;
+             break;
            }
          }
          if (!digitExists) {
-           AliMFTDigit newSDigit;
-           newSDigit.SetEloss(0.);
-           newSDigit.SetDetElemID(hit->GetDetElemID());
-           newSDigit.SetPlane(hit->GetPlane());
-           newSDigit.AddMCLabel(hit->GetTrack());
-           newSDigit.SetPixID(xPixelNew, yPixelNew, 0);
-           newSDigit.SetPixWidth(fSegmentation->GetPixelSizeX(sDigit.GetDetElemID()), 
-                                 fSegmentation->GetPixelSizeY(sDigit.GetDetElemID()),
-                                 fSegmentation->GetPixelSizeZ(sDigit.GetDetElemID()));  
-           newSDigit.SetPixCenter(fSegmentation->GetPixelCenterX(sDigit.GetDetElemID(), xPixelNew), 
-                                  fSegmentation->GetPixelCenterY(sDigit.GetDetElemID(), yPixelNew),
-                                  fSegmentation->GetPixelCenterZ(sDigit.GetDetElemID(), 0)); 
-           new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(newSDigit);
+           AliMFTDigit *sideSDigit = new AliMFTDigit();
+           sideSDigit->SetEloss(0.);
+           sideSDigit->SetDetElemID(hit->GetDetElemID());
+           sideSDigit->SetPlane(hit->GetPlane());
+           sideSDigit->AddMCLabel(hit->GetTrack());
+           sideSDigit->SetPixID(xPixel, yPixel, 0);
+           sideSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(sideSDigit->GetDetElemID()), 
+                                   fSegmentation->GetPixelSizeY(sideSDigit->GetDetElemID()),
+                                   fSegmentation->GetPixelSizeZ(sideSDigit->GetDetElemID()));  
+           sideSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(sideSDigit->GetDetElemID(), xPixel), 
+                                    fSegmentation->GetPixelCenterY(sideSDigit->GetDetElemID(), yPixel),
+                                    fSegmentation->GetPixelCenterZ(sideSDigit->GetDetElemID(), 0)); 
+           new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*sideSDigit);
+         }
+       }
+      }
+    }
+    
+    // ------------ In case we should simulate a rectangular pattern of pixel...
+    
+    if (fSegmentation->GetPlane(mainSDigit->GetPlane())->HasPixelRectangularPatternAlongY()) {
+      for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
+       AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
+       if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {   // both pair or both odd
+         xPixel = mySDig->GetPixelX();
+         yPixel = mySDig->GetPixelY()+1;
+         if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
+           AliMFTDigit *newSDigit = new AliMFTDigit();
+           newSDigit->SetEloss(0.);
+           newSDigit->SetDetElemID(mySDig->GetDetElemID());
+           newSDigit->SetPlane(mySDig->GetDetElemID());
+           newSDigit->SetPixID(xPixel, yPixel, 0);
+           newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
+                                  fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
+                                  fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
+           newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
+                                   fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
+                                   fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0)); 
+           new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
+         }
+       }
+       else {   // pair-odd
+         xPixel = mySDig->GetPixelX();
+         yPixel = mySDig->GetPixelY()-1;
+         if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
+           AliMFTDigit *newSDigit = new AliMFTDigit();
+           newSDigit->SetEloss(0.);
+           newSDigit->SetDetElemID(mySDig->GetDetElemID());
+           newSDigit->SetPlane(mySDig->GetPlane());
+           newSDigit->SetPixID(xPixel, yPixel, 0);
+           newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
+                                  fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
+                                  fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
+           newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
+                                   fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
+                                   fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0)); 
+           new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
          }
        }
       }
     }
 
-    for (Int_t iSideDigit=0; iSideDigit<fSideDigits->GetEntries(); iSideDigit++) {
-      AliMFTDigit *newDigit = (AliMFTDigit*) fSideDigits->At(iSideDigit);
-      new ((*pSDigList[sDigit.GetPlane()])[pSDigList[sDigit.GetPlane()]->GetEntries()]) AliMFTDigit(*newDigit);
+    // -------- checking which pixels switched on have their diode actually within the charge dispersion radius
+
+    for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
+      AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
+      Double_t distance = TMath::Sqrt(TMath::Power(mySDig->GetPixelCenterX()-hit->X(),2) + TMath::Power(mySDig->GetPixelCenterY()-hit->Y(),2));
+      if (fSegmentation->GetPlane(mySDig->GetPlane())->HasPixelRectangularPatternAlongY()) {
+       if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {  // both pair or both odd
+         if (distance<fChargeDispersion) {
+           AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
+                            mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
+           new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
+           xPixel = mySDig->GetPixelX();
+           yPixel = mySDig->GetPixelY()+1;
+           if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
+             AliMFTDigit *newSDigit = new AliMFTDigit();
+             newSDigit->SetEloss(0.);
+             newSDigit->SetDetElemID(mySDig->GetDetElemID());
+             newSDigit->SetPlane(mySDig->GetPlane());
+             newSDigit->SetPixID(xPixel, yPixel, 0);
+             newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
+                                    fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
+                                    fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
+             newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
+                                     fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
+                                     fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0));
+             AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
+                              newSDigit->GetPixelCenterX(), newSDigit->GetPixelCenterY(), newSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
+             new ((*pSDigList[newSDigit->GetPlane()])[pSDigList[newSDigit->GetPlane()]->GetEntries()]) AliMFTDigit(*newSDigit);
+           }
+         }
+       }
+      }
+      else {
+       if (distance<fChargeDispersion) {
+         AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
+                          mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
+         new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
+       }
+      }
     }
 
-    fSideDigits->Clear();    
+    fSideDigits->Delete(); 
 
   }
-  
-  AliDebug(1,"Stop Hits2SDigitsLocal");
+
+  AliDebug(1,"Exiting Hits2SDigitsLocal");
 
 }
 
@@ -507,8 +605,6 @@ void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t
 
 void AliMFT::MakeBranch(Option_t *option) {
 
-  printf("AliMFT::MakeBranch(...)\n");
-
   // Create Tree branches 
   AliDebug(1, Form("Start with option= %s.",option));
   
@@ -590,7 +686,7 @@ void AliMFT::SetTreeAddress() {
 
 void AliMFT::SetGeometry() {
 
-  printf("AliMFT::SetGeometry\n");
+  AliInfo("AliMFT::SetGeometry\n");
 
   fSegmentation = new AliMFTSegmentation(fNameGeomFile.Data());