Logging of Debug, Info and Error Messages follwing AliRoot Standard http://aliweb...
[u/mrichter/AliRoot.git] / MUON / AliMUONproto.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 /*
19 Implementation of MUON Chamber Prototype Class 
20 Both read digits from raw data or use the Monte-Carlo. 
21 1-Feb-2000 Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
22 */
23
24
25 ////////////////////////////////////////////////
26 //  Manager and hits classes for set:PROTO    //
27 ////////////////////////////////////////////////
28
29 #include <TTUBE.h>
30 #include <TBRIK.h>
31 #include <TRotMatrix.h>
32 #include <TNode.h> 
33 #include <TTree.h> 
34 #include <TRandom.h> 
35 #include <TObject.h>
36 #include <TVector.h>
37 #include <TObjArray.h>
38 #include <TMinuit.h>
39 #include <TParticle.h>
40 #include <TROOT.h>
41 #include <TFile.h>
42 #include <TNtuple.h>
43 #include <TCanvas.h>
44 #include <TPad.h>
45 #include <TDirectory.h>
46 #include <TObjectTable.h>
47 #include <AliPDG.h>
48 #include <TArrayI.h>
49 #include <AliDetector.h>
50
51 #include "AliMUONChamber.h"
52 #include "AliMUONproto.h"
53 #include "AliMUONHit.h"
54 #include "TTUBE.h"
55 #include "AliMUONClusterFinder.h"
56 #include "AliRun.h"
57 #include "Riostream.h"
58 #include "AliCallf77.h" 
59 #include "AliConst.h"
60 //#include "chainalice2.h"
61 #include "AliMUONSegmentationV0.h"
62 //#include "AliMUONSegResV11.h"
63 #include "AliLog.h"
64
65 ClassImp(AliMUONproto)
66
67 //___________________________________________    
68 AliMUONproto::AliMUONproto()
69         : AliMUON()
70 {
71     AliInfo(" Calling AliMUONproto constructor...");
72     
73         //
74         //
75         //
76 }
77  
78 //___________________________________________
79 AliMUONproto::AliMUONproto(const char *name, const char *title) 
80         : AliMUON(name,title)
81 //  AliMUON defines everything, but the chamber for NCH=1
82 {
83     
84 //  z-Positions of Chambers
85     const Float_t zch[1] = {975.};
86 //
87 //  inner diameter
88     const Float_t dmi[1] = {0.};
89 //
90 //  outer diameter
91     const Float_t dma[1] = {29.2};
92 //
93 //    
94 //    Default Parameters for ALICE2 prototype
95
96 //
97     (*fChambers)[0] = new AliMUONChamber();
98     AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[0];
99     chamber->SetGid(0);
100     chamber->SetZ(zch[0]);
101 //
102     chamber->InitGeo(zch[0]); 
103     chamber->SetRInner(dmi[0]/2);
104     chamber->SetROuter(dma[0]/2);
105
106     for (int i = 0; i <= 99; i++) {
107         fThreshold[i] = 0.;
108     }
109     
110 }
111
112 #ifdef WE_FORGET_FOR_THE_MOMENT
113 //___________________________________________ 
114 void AliMUONproto::GetRawDigits(Int_t evnb, Int_t *lptr, Int_t ilen) {
115
116     Int_t ip = 0;
117     Int_t equip = 0;
118     Int_t nochnl;
119     Int_t loop;
120     Int_t val;    
121     Int_t itype;
122     Int_t id;
123     Int_t serial;
124     Int_t equiplen;
125
126     Int_t digits[5]; 
127     
128     AliMUON* MUON  = (AliMUON*)gAlice->GetModule("MUON");
129     AliMUONSegmentationV0* seg = (AliMUONSegmentationV0*) Chamber(0).SegmentationModel(1);
130     
131     Int_t tracks[10];
132     Int_t charges[10];
133         
134     for (Int_t i = 0; i < 10; i++) {
135         tracks[i] = 0;
136         charges[i] = 0;
137     }
138         
139     Int_t ich = 0;
140     
141   nwtype:             
142
143     itype = lptr[ip++];
144     id = lptr[ip++];
145     equiplen = lptr[ip++];
146     
147     if (equiplen < 0 ) return;
148
149     if (itype != (int)(0XCACCA008)) {        
150         ip += equiplen;
151         if (ip < ilen-2)  goto nwtype;
152     }    
153     else {
154         serial = id >> 16;
155         equip = id & 0x1;
156         if ((serial == 190) && (equip == 1)) {
157             for (loop = 0; loop < equiplen; loop++) {
158                 nochnl = (lptr[ip] & 0x7ff000 ) >> 12;
159                 val = lptr[ip] & 0x3ff;
160                 // fill digits from raw data according to cathode connexions
161                 if (group[nochnl][2][1]!=0)
162                     digits[0] = group[nochnl][2][1] - 12;
163                 else if  (group[nochnl][1][1]!=0)
164                     digits[0] = group[nochnl][1][1] - 12;
165                 else
166                     digits[0] = group[nochnl][0][1] - 12;
167                 digits[1] = group[nochnl][0][0];
168                 if (digits[0] != seg->Ix(digits[0], digits[1]))
169                     printf("Pb pour ix=%d,iy=%d\n", digits[0], digits[1]);
170                 digits[2] = val;
171                 digits[3] = 0;
172                 digits[4] = 0;
173                 if (digits[2] >= fThreshold[nochnl]) 
174                     MUON->AddDigits(ich, tracks, charges, digits);
175
176                 ip++;
177             }   
178         } 
179         else
180             ip += equiplen;
181         
182         if (ip < ilen-2) goto nwtype;
183     }
184     
185     gAlice->TreeD()->Fill();
186     MUON->ResetDigits();
187
188     gAlice->TreeD()->Fill();
189     MUON->ResetDigits();
190     
191     char hname[30];
192     sprintf(hname, "TreeD%d", evnb);
193     gAlice->TreeD()->Write(hname,TObject::kOverwrite);
194     // reset tree
195     gAlice->TreeD()->Reset();
196
197 }
198 #endif
199
200 //___________________________________________
201 void AliMUONproto::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
202 {
203     Int_t i=2*(id-1);
204     ((AliMUONChamber*)(*fChambers)[i])->SetPadSize(isec,p1,p2);
205 }
206
207 //___________________________________________
208 void AliMUONproto::SetChargeSlope(Int_t id, Float_t p1)
209 {
210     Int_t i=2*(id-1);
211     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
212 }
213
214 //___________________________________________
215 void AliMUONproto::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
216 {
217     Int_t i=2*(id-1);
218     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
219 }
220
221 //___________________________________________
222 void AliMUONproto::SetMaxAdc(Int_t id, Float_t p1)
223 {
224     Int_t i=2*(id-1);
225     ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
226 }
227
228 //___________________________________________
229 void AliMUONproto::CreateGeometry()
230 {
231     Int_t *idtmed = fIdtmed->GetArray()-1099;
232 //
233     Float_t zpos;
234     Float_t bpar[3];
235     Int_t idGas=idtmed[1100];
236     
237     AliMUONChamber *iChamber;
238 //********************************************************************
239 //                          Prototype ALICE2                        **
240 //********************************************************************
241     iChamber=(AliMUONChamber*) (*fChambers)[0];
242     
243     zpos=iChamber->Z(); 
244 //
245     const Float_t X_POS = 11*.975/2; //half x
246     const Float_t Y_POS = 18*.55/2;  //half y
247     const Float_t Z_POS = 0.325;
248     
249     bpar[0] = X_POS; 
250     bpar[1] = Y_POS;  
251     bpar[2] = Z_POS;
252     
253     gMC->Gsvolu("C01G", "BOX", idGas, bpar, 3);
254     gMC->Gspos("C01G", 1, "ALIC", -bpar[0], bpar[1], zpos, 0, "ONLY");
255 }
256
257 //___________________________________________
258 void AliMUONproto::CreateMaterials()
259 {
260   // *** DEFINITION OF AVAILABLE MUON MATERIALS *** 
261   //
262   //     Ar-CO2 gas 
263     Float_t ag1[3]   = { 39.95,12.01,16. };
264     Float_t zg1[3]   = { 18.,6.,8. };
265     Float_t wg1[3]   = { .8,.0667,.13333 };
266     Float_t dg1      = .001821;
267     
268     Float_t epsil  = .001; // Tracking precision, 
269     Float_t tmaxfd = -20.; // Maximum angle due to field deflection
270     Float_t stmin  = -.8;
271     Int_t ISXFLD   = gAlice->Field()->Integ();
272     Float_t SXMGMX = gAlice->Field()->Max();
273
274     //
275     // --- Define the various materials for GEANT ---
276     AliMixture(22, "ArCO2 80%$", ag1, zg1, dg1, 3, wg1);
277
278     //
279     //    Ar/CO2
280     AliMedium(1, "ARG_CO2   ", 22, 1, ISXFLD, SXMGMX, tmaxfd, fMaxStepGas,
281               fMaxDestepAlu, epsil, stmin);
282     //    Air 
283         //AliMedium(1, "AIR_CH_US         ", 15, 1, ISXFLD, SXMGMX, tmaxfd, -1., -.3, epsil, stmin);
284 }
285
286 //___________________________________________
287 void AliMUONproto::Init()
288 {
289    printf("\n\n\n Start Init for Prototype ALICE2 - CPC chamber type\n\n\n");
290
291    // 
292    // Initialize Tracking Chambers
293    //
294    
295    ( (AliMUONChamber*) (*fChambers)[0])->Init();
296    
297    //
298    // Set the chamber (sensitive region) GEANT identifier
299    ((AliMUONChamber*)(*fChambers)[0])->SetGid(gMC->VolId("C01G"));
300
301    printf("\n\n\n Finished Init for Prototype ALICE2 - CPC chamber type\n\n\n");
302
303 }
304
305 //___________________________________________
306 void AliMUONproto::StepManager()
307 {
308   Int_t          copy, id;
309   static Int_t   idvol;
310   static Int_t   vol[2];
311   Int_t          ipart;
312   TLorentzVector pos;
313   TLorentzVector mom;
314   Float_t        theta,phi;
315   Float_t        destep, step;
316   
317   static Float_t eloss, eloss2, xhit, yhit, tlength;
318   const  Float_t big=1.e10;
319   
320   //  modifs perso
321   static Float_t hits[14];
322
323   TClonesArray &lhits = *fHits;
324
325   //
326   // Set maximum step size for gas
327   // numed=gMC->GetMedium();
328   //
329   // Only charged tracks
330   if( !(gMC->TrackCharge()) ) return; 
331   //
332   // Only gas gap inside chamber
333   // Tag chambers and record hits when track enters 
334   idvol=-1;
335   id=gMC->CurrentVolID(copy);
336   
337     for (Int_t i=1; i<=kNCH; i++) {
338       if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()){ 
339           vol[0]=i; 
340           idvol=i-1;
341       }
342     }
343     if (idvol == -1) return;
344   //
345   // Get current particle id (ipart), track position (pos)  and momentum (mom) 
346   gMC->TrackPosition(pos);
347   gMC->TrackMomentum(mom);
348
349   ipart  = gMC->TrackPid();
350   //Int_t ipart1 = gMC->IdFromPDG(ipart);
351   //printf("ich, ipart %d %d \n",vol[0],ipart1);
352
353   //
354   // momentum loss and steplength in last step
355   destep = gMC->Edep();
356   step   = gMC->TrackStep();
357   
358   //
359   // record hits when track enters ...
360   if( gMC->IsTrackEntering()) {
361       gMC->SetMaxStep(fMaxStepGas);
362       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
363       Double_t rt = TMath::Sqrt(tc);
364       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
365       Double_t tx=mom[0]/pmom;
366       Double_t ty=mom[1]/pmom;
367       Double_t tz=mom[2]/pmom;
368       Double_t s=((AliMUONChamber*)(*fChambers)[idvol])
369           ->ResponseModel()
370           ->Pitch()/tz;
371       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
372       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
373       hits[0] = Float_t(ipart);         // Geant3 particle type
374       hits[1] = pos[0]+s*tx;                 // X-position for hit
375       hits[2] = pos[1]+s*ty;                 // Y-position for hit
376       hits[3] = pos[2]+s*tz;                 // Z-position for hit
377       hits[4] = theta;                  // theta angle of incidence
378       hits[5] = phi;                    // phi angle of incidence 
379       hits[8] = (Float_t) fNPadHits;   // first padhit
380       hits[9] = -1;                     // last pad hit
381
382       // modifs perso
383       hits[10] = mom[3]; // hit momentum P
384       hits[11] = mom[0]; // Px/P
385       hits[12] = mom[1]; // Py/P
386       hits[13] = mom[2]; // Pz/P
387       // fin modifs perso
388
389       // phi angle of incidence
390       tlength = 0;
391       eloss   = 0;
392       eloss2  = 0;
393       xhit    = pos[0];
394       yhit    = pos[1];      
395       // Only if not trigger chamber
396       if(idvol<10) {
397           //
398           //  Initialize hit position (cursor) in the segmentation model 
399           ((AliMUONChamber*) (*fChambers)[idvol])
400               ->SigGenInit(pos[0], pos[1], pos[2]);
401       } else {
402           //geant3->Gpcxyz();
403           //printf("In the Trigger Chamber #%d\n",idvol-9);
404       }
405   }
406   eloss2+=destep;
407   
408   // 
409   // Calculate the charge induced on a pad (disintegration) in case 
410   //
411   // Mip left chamber ...
412   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
413       gMC->SetMaxStep(big);
414       eloss   += destep;
415       tlength += step;
416       
417       // Only if not trigger chamber
418       if(idvol<10) {
419           if (eloss > 0) 
420               MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
421       }
422           
423       hits[6]=tlength;
424       hits[7]=eloss2;
425       if (fNPadHits > (Int_t)hits[8]) {
426           hits[8]= hits[8]+1;
427           hits[9]= (Float_t) fNPadHits;
428       }
429     
430       new(lhits[fNhits++]) 
431           AliMUONHit(fIshunt,gAlice->GetCurrentTrackNumber(),vol,hits);
432       eloss = 0; 
433       //
434       // Check additional signal generation conditions 
435       // defined by the segmentation
436       // model (boundary crossing conditions) 
437   } else if 
438       (((AliMUONChamber*) (*fChambers)[idvol])
439        ->SigGenCond(pos[0], pos[1], pos[2]))
440   {
441       ((AliMUONChamber*) (*fChambers)[idvol])
442           ->SigGenInit(pos[0], pos[1], pos[2]);
443 //      printf("\n-> MakePadHits, reason special %d",ipart);
444       if (eloss > 0)
445           MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
446       xhit     = pos[0];
447       yhit     = pos[1]; 
448       eloss    = destep;
449       tlength += step ;
450       //
451       // nothing special  happened, add up energy loss
452   } else {        
453       eloss   += destep;
454       tlength += step ;
455   }
456 }
457
458 //___________________________________________
459 void AliMUONproto::BuildGeometry()
460 {
461     TNode *Node;
462     TNode *Top;
463     
464     const int kColorMUON = kBlue;
465     //
466     Top=gAlice->GetGeometry()->GetNode("alice");
467     //
468     //
469     //
470     float dx, dy, dz, zpos;
471
472     const Float_t cz[1]={975.};
473     
474     zpos=cz[0];
475
476     dx=11*.975/2;
477     dy=18*.55/2;
478     dz=0.325;
479
480     new TBRIK("C_MUON101","Mother volume for Proto.","void",dx*2,dy*2,dz*2);
481     TBRIK* PROTO = new TBRIK("PROT", "Proto. Sens. Region","void",dx,dy,dz);
482     Top->cd();
483     Node = new TNode("MUON101","ChamberNode","C_MUON101",0,0,zpos,"");
484     Node->SetLineColor(kColorMUON);
485     Node->SetVisibility(0);    
486     fNodes->Add(Node);
487     Node->cd();
488     Node = new TNode("MUON201", "Proto. Sens. Region Node", PROTO, -dx, dy, dz);
489     Node->SetLineColor(kColorMUON);    
490
491 }
492
493 //_____________________________________________________________________________
494 void AliMUONproto::FindClusters(Int_t nev, Int_t last_entry)
495 {
496
497 //
498 // Loop on chambers and on cathode planes
499 //
500   for (Int_t icat = 0; icat < 2; icat++) {
501             gAlice->ResetDigits();
502             gAlice->TreeD()->GetEvent(last_entry+icat);
503             
504       for (Int_t ich = 0; ich < kNTrackingCh; ich++) {
505           AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
506           TClonesArray *MUONdigits  = this->DigitsAddress(ich);
507           if (MUONdigits == 0) continue;
508           //
509           // Get ready the current chamber stuff
510           //
511
512           AliMUONResponse* response = iChamber->ResponseModel();
513           AliMUONSegmentation* seg = iChamber->SegmentationModel(icat+1);
514           AliMUONClusterFinder* rec = iChamber->ReconstructionModel();
515
516           if (seg) {      
517               rec->SetSegmentation(seg);
518               rec->SetResponse(response);
519               rec->SetDigits(MUONdigits);
520               rec->SetChamber(ich);
521               rec->FindRawClusters();
522           }
523           
524           TClonesArray *fRch;
525           fRch=RawClustAddress(ich);
526           fRch->Sort();
527           // it seems to work
528           
529       } // for ich
530       // fill the tree
531       gAlice->TreeR()->Fill();
532
533       ResetRawClusters();
534   } // for icat
535
536   char hname[30];
537   sprintf(hname,"TreeR%d",nev);
538   gAlice->TreeR()->Write(hname);
539   gAlice->TreeR()->Reset();
540 }
541
542 #ifdef WE_FORGRET_THIS_SHIT
543 //_____________________________________________________________________________
544 void AliMUONproto::SetThreshold()
545 {
546
547     ifstream inputFile("/home/alice/guernane/aliroot/pro/MUON/crped190.dat", ios::in);
548
549     if (inputFile.fail()) {
550         AliError("Error opening file");
551         exit(2);
552     }
553
554     char buff[32];
555     Int_t Serial;
556     Int_t Ntrigger;
557     Int_t Nchannel;
558     Int_t i1;
559     Int_t i2;
560         
561     inputFile >> buff;
562
563     inputFile >> Serial;
564     inputFile >> Ntrigger;
565     inputFile >> Nchannel;
566     inputFile >> i1;
567     inputFile >> i2;
568
569     Float_t ped0[Nchannel];
570     Float_t sig0[Nchannel];
571     Float_t ped1[Nchannel];
572     Float_t sig1[Nchannel];
573     Int_t ichannel;
574
575     for (Int_t i = 0; i < Nchannel-1; i++) {
576     ped0[i] = 0;
577     sig0[i] = 0;
578     ped1[i] = 0;
579     sig1[i] = 0;        
580     }
581     
582     for (Int_t i = 0; i < Nchannel-1; i++) {
583         inputFile >> ichannel;
584         inputFile >> ped0[i];
585         inputFile >> sig0[i];
586         inputFile >> ped1[i];
587         inputFile >> sig1[i];
588         fThreshold[i] = fNsigma*sig1[i];        
589     }
590     
591     inputFile.close();    
592 }
593 #endif